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ABSTRACT 


Requirements  for  submarine  periscope  depth  operations  have  been  increased  by 
integration  with  carrier  battle  groups,  littoral  operations,  and  contributions  to  joint  surveillance. 
Improved  periscope  depth  performance  is  therefore  imperative.  Submarine  control  personnel 
rely  on  a  large  number  of  analog  gauges  and  indications.  An  integrated  digital  display  system 
could  enhance  the  ergonomics  of  the  human  control  interface  and  display  additional 
parameters.  This  thesis  investigates  the  required  feedbacks  for  robust  automatic  depth  control 
at  periscope  depth,  and  thus  indirectly  determines  the  additional  parameters  desired  for  an 
integrated  display. 

A  model  of  vertical  plane  submarine  dynamics  is  coupled  with  first  and  second  order 
wave  force  solutions  for  a  particular  submarine  hull  form.  Sliding  mode  control  and  several 
schemes  of  state  feedback  are  used  for  automatic  control.  Head  and  beam  seas  at  sea  states 
three  and  four  are  investigated.  The  automatic  control  effectiveness  provides  insight  into  the 
indications  used  by  the  ship's  control  party  in  operations  at  periscope  depth.  One  possible 
display  system  is  proposed,  with  several  additional  enhancements  to  improve  ship's  safety, 
reduce  operator  fatigue,  and  enable  accurate  reconstmction  of  the  events  leading  to  a  loss  of 
depth  control. 
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1.  INTRODUCTION 


A,  GENERAL 

The  need  for  attack  submarines  to  operate  at  periscope  depth  has  been  increased  by 
integration  with  carrier  battle  groups,  operations  in  the  shallow  littoral,  and  contributions  to 
joint  surveillance. 

Operating  at  periscope  depth  beneath  a  seaway,  a  submarine  is  in  an  unstable 
condition.  As  the  free  surface  is  approached,  the  seaway  forces  increase,  trying  to  pull  the 
submarine  to  the  surface.  To  counter  these  forces,  the  ship’s  ballast  is  changed  and  control 
surfaces  are  used.  Because  of  the  seaways  stochastic  nature,  manual  operation  for  long  periods 
at  periscope  depth  taxes  the  ship’s  control  party. 

Operators  must  remain  aware  of  the  environmental  conditions.  If  the  sea  becomes 
quiescent,  the  submarine  will  sink  out.  If  the  sea  suction  forces  are  greater  than  the  ballast  and 
planes  authority,  the  submarine  will  broach  the  free  surface  increasing  detection  risk  by  several 
orders  of  magnitude.  Other  events,  like  temperature  or  salinity  changes,  can  also  have  major 
effects  on  reliable  depth  keeping.  Contributing  to  the  environmental  issues,  the  need  to  use 
minimum  speed  for  a  given  sea  state  to  control  the  detectable  mast  feather  reduces  the 
available  planes  authority,  and  increases  the  difficulty  of  depth  control. 

However,  the  current  submarine  force  is  not  optimized  for  these  operations.  One 
inexpensive  area  for  improvement  is  the  display  system  for  the  ship’s  control  party.  Modern 
digital  display  systems  offer  ergonomic  improvements  over  current  gauges  and  readouts. 

Given  a  requirement  to  conduct  submarine  ship  control  manually,  a  fundamental 
question  is  that  of  how  to  display  the  state  of  the  ship  to  the  operators.  Aside  from  the 
obvious  indications  like  ship’s  pitch  angle,  depth,  and  control  surface  positions  what  else  would 
be  useful?  Candidates  include  the  net  force  acting  on  the  ship,  accelerations,  and  various  time 
averaged  values.  Implied  in  this  is  that  a  nontraditional  means  of  display  will  be  used  to  show 
these  parameters,  so  that  the  operators  will  not  have  to  rely  on  a  number  of  gauges  or  meters, 
with  averaging  of  results  only  available  only  by  the  calibrated  eye. 

An  intelligent  assistant  to  the  ship’s  control  party  would  show  items  of  current 
concern,  and  issue  alerts  based  on  an  operator  programmable  doctrine.  Issues  like  mast 
exposure,  ship’s  relationship  to  the  bottom,  and  trim  state  could  be  shown  in  an  intuitive, 
logical  manner. 
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Current  evolutions  and  other  items  relating  to  the  tactical  employment  could  be  included  as 
required. 


B.  AIM  OF  THIS  STUDY 

Although  the  ship’s  control  party  currendy  relies  on  a  small  number  of  indications,  the 
ability  to  sense  “by  the  seat  of  the  pants”  cannot  be  discounted.  This  thesis  investigates 
required  feedbacks  for  robust  automatic  depth  control  at  periscope  depth,  and  thus  indirecdy 
evaluates  the  additional  indications  to  be  added  to  an  integrated  display. 

This  approach  assumes  that  the  best  ship’s  control  parties  already  use  system  states  for 

control  which  are  not  explicidy  displayed. 

C.  THESIS  OUTLINE 

Chapter  II  contains  the  development  of  the  deeply  submerged  submarine  dynamics 
model.  Chapter  III  gives  the  development  and  source  of  the  wave  forces  used  to  simulate 
operations  at  periscope  depth.  In  Chapter  IV,  optimization  studies  are  performed  for  nine 
different  cases  of  state  feedback  control.  This  gives  a  feeling  for  the  quality  of  depth  control 
achievable  by  the  use  of  different  levels  of  sensors.  Chapter  V  explores  the  use  of  sliding 
mode  control  for  periscope  depth  operations.  In  Chapter  VI,  current  ships  control  technology 
is  reviewed  and  an  integrated  display  is  proposed.  Conclusions  and  recommendations  are 
given  in  Chapter  VII. 


2 


II.  SUBMARINE  DYNAMICS  MODEL 


A.  INTRODUCTION 

When  a  submarine  is  deeply  submerged,  many  of  its  maneuvering  characteristics  can 
be  determined  from  application  of  Morison’s  equation  to  model  test  data.  A  series  of  trials, 
often  done  with  a  planar  motion  mechanism  (PMM),  give  the  damping  and  inertia  coefficients 
for  small  maneuvers  in  each  of  the  six  degrees  of  freedom.  This  method  is  not  without  limits. 
For  trials  done  in  the  horizontal  and  vertical  planes  only,  nonlinear  cross  coupling  effects  are 
ignored.  The  hydrodynamic  coefficients  work  poorly  for  prediction  of  high  speed  maneuvers 
and  control  surface  casualties.  Here  the  large  crossflow  velocities,  vortex  hull  interaction,  and 
flow  separation  all  have  effects  which  are  not  predicted  by  the  hydrodynamic  coefficients.  It  is 
possible,  however,  to  include  some  of  these  effects  as  additional  nonlinear  terms. 

As  the  submarine  approaches  the  free  surface,  several  complexities  are  introduced  into 
the  hydrodynamic  coefficient  approach.  First,  the  inertia  terms  change  as  an  acceleration  will 
no  longer  act  upon  an  effectively  infinite  region.  Second,  an  inviscid  form  of  damping  exists 
near  the  free  surface.  This  comes  about  from  the  generation  of  waves  by  the  body,  and 
depends  on  the  body  depth  and  character  of  motion.  Finally,  the  interaction  between  the 
incident  waves  and  the  submarine  introduces  added  forces  and  moments.  These  effects 
combine  to  make  designing  for  periscope  depth  vexing  for  engineers  and  operating  at 
periscope  depth  an  art  for  the  ship’s  crew. 

The  approach  in  this  thesis  will  be  to  first  establish  a  dynamics  model  appropriate  for  a 
deeply  submerged  submarine  at  low  to  moderate  speeds.  The  forces  and  moments  resulting 
from  the  seaway  will  then  be  superimposed  on  this  model  to  provide  a  reasonable 
approximation  to  the  submarine  motion  beneath  waves. 

B.  DEEPLY  SUBMERGED  EQUATIONS  OF  MOTION 

1.  Definition  of  coordinate  system  and  states 

The  coordinate  system  defined  in  Figure  1  will  be  used.  The  origin  of  the  global 
coordinate  system  is  fixed  at  the  ocean  surface.  The  ^  axis  is  positive  downward,  towards  the 
ocean  bottom.  The  x  axis  is  positive  in  the  direction  of  intended  submarine  motion.  The  body 
fixed  coordinates  are  rotated  from  the  global  coordinates  by  the  angle  6 .  Body  fixed  velocities 
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w  (heave)  ,  u  (surge),  and  q  (pitch)  are  shown.  The  control  surface  deflections,  Sj,  (bow  planes) 
and  5^  (stern  planes)  are  also  defined. 


2.  Hydrodynamic  coefficients  review 

For  a  deeply  submerged  submarine,  small  motions  can  be  analyzed  using  the  concept 
of  hydrodynamic  coefficients.  These  represent  a  Taylor  series  expansion  of  the  functional 
relationship  between  body  movements  and  the  resulting  fluid  forces.  For  example,  given  the 
deeply  submerged  body  in  Figure  2  undergoing  pure  heave,  resulting  body  forces  can  be 
expressed  in  the  following  manner: 
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Figure  2.  Submerged  body  in  pure  heave 

Af  =  M^w+ W 

Z  =  Zy,M>+Z^^w\\^  +  Z^w  (2) 

This  method  is  extended  to  the  six  degrees  of  freedom  of  the  body,  and  done  for 
velocity  and  acceleration  components  of  the  movement.  This  includes  representations  of 
added  mass,  viscous  drag,  and  square  law  drag. 

3.  Vertical  plane  equations  of  motion 

By  using  this  system  of  notation,  and  applying  Newton’s  second  law  to  the  body  fixed 
coordinates,  and  transforming  to  global  coordinates,  the  equations  of  pitch  and  heave  may  be 
obtained  in  the  vertical  plane.  The  general  case  is  quite  complex,  having  centers  of  mass  and 
buoyancy  that  are  separate  firom  each  other  and  the  coordinate  system  origin.  This,  along  with 
cross  coupled  hydrodynamic  coefficients,  results  in  a  nonlinear,  coupled  set  of  differential 
equations. 

These  equations  of  pitch  and  heave  may  be  simplified  considerably  by  several 
reasonable  assumptions.  Assuming  that  the  submarine  motion  is  constrained  to  the  vertical 
plane,  the  equations  of  motion  for  heave  and  pitch  are  (Smith,  Crane,  and  Summey  (1978)): 

rriw-uq-XQq-Zcq^]=  P) 

Z^q  +  Z^,w  +  Z^uq  +  Z„w 

+  u\Zs5,+Zs5,) 
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(4) 


Iyq-nixc{w-uq)-Za(.u  +  '<^^)]=^ 

M^q+M^w+  M^uq+M^uw 

+  u^(Mg^Si,  +  Mg5^) 

-  (x^mg  -  Xg  8)005(6) 

-(Zaing- Zb  B)sin(d) 

It  is  apparent  that  Equations  (3)  and  (4)  are  nonlinear,  coupled  differential  equations  in 
u>  and  q  and  ».  To  reduce  this  coupling,  terms  involving  the  derivatives  of  o'  and  q  can  be 
collected,  resulting  in  a  mass  matrix. 

_  m-Z^  -Z^  (5) 

The  mass  matrix  can  be  readily  inverted: 

z,  1  (6) 

_  ,  "'-Zw 

= - = - - — 

{m-Z^)(Iy  -  M^)-ZgM^ 

By  applying  Equation  (6),  the  cross  coupling  of  terms  in  vv  and  q  can  be  removed 
from  Equations  (3)  and  (4).  To  allow  the  introduction  of  external  forces  and  moments,  the 
system  was  augmented  by  force  and  moment  disturbances  acting  at  the  origin  of  the  body 
fixed  coordinates.  They  were  multiplied  by  the  cosine  of  the  pitch  angle  for  conversion  to  the 
body  fixed  coordinate  system.  These  disturbances  can  be  used  to  input  external  effects,  such 
as  changes  in  trim  and  wave  forces.  By  further  assuming  that  the  center  of  buoyancy  is  at  the 
body  fixed  coordinate  system  origin,  the  center  of  mass  is  directly  below,  and  that  the  forward 
speed  u  is  constant,  the  equations  of  motion  can  be  reduced  to  the  following: 

\v  =  a]]Uw  +  a]2uq+a^■y5in(6)+b]]U^5f,  +fc,2M^5,  +  cos(0)  +  ei]^^  +e,2?>v 

q  =  a2]Uw  +  a22uq  +U23  sin(0)  +^>21“^^*  ^t/  cos(0)  +  Oi\q  +  ^22^^ 

e  =  q 

z  =  wcos(0)-«sin(0) 

X  -  wsini^)  +  Mcos(0) 


(7) 

(8) 
(9) 

(10) 

(11) 
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where: 


_  iZ,+m)U,-  M-)  +  Z.M^ 
(m-Z^.)(/,  - 

M^Zi,  +  (m-Z^)M„ 
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Z,)M^ 
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_ Z^Zgin 
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{m-ZJZr,m 

Equations  (7)  through  (1 1)  are  the  governing  equations  of  motion  for  this  thesis.  It  is 
of  note  that  the  disturbance  force  and  moment  terms  represent  accelerations  due  to  the 
disturbances.  To  provide  ease  of  use,  the  equations  of  motion  were  implemented  in  the 
SIMULINK®  model  shown  in  Figure  3.  This  building  block  approach  was  very  effective  for 
conducting  studies  on  the  effectiveness  of  different  types  of  controllers. 


Figure  3.  SIMULINK®  model  of  vertical  plane  submarine  dynamics 

For  control  design,  it  is  convenient  to  use  a  linear  state  space  representation  of  the 
system.  This  allows  the  use  of  a  variety  of  controller  design  tools  including  pole  placement  and 
linear  quadratic  regulator  algorithms.  Equations  (7)  through  (11)  can  be  linearized  about  a  level 
flight  condition.  This  results  in  the  linear  state  space  representation: 


w  =  aiiuw  a^2^q  +  a^-^6  + 

(12) 

q  =  a2\UW  +  a22^^  ^ 

(13) 

e  =  q 

(14) 

Z^w-uO 

(15) 

X  =  w6  +  u 

(16) 
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Equations  (12)  through  (15)  can  be  rewritten  in  matrix  form.  This  form  of  the  linear 
submarine  vertical  plane  dynamics  equations  will  be  used  for  controller  design.  For  controller 
design,  Equation  (16)  was  excluded  from  the  matrix  form.  Because  of  the  constant  forward 
speed  u  assumption,  there  was  no  direct  means  of  control  for  x. 
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C.  EXTENSION  TO  VERTICAL  PLANE  PATHKEEPING 

Equations  (7)  through  (10)  and  the  corresponding  SIMUUNK®  model  are  linearized 
around  a  constant  commanded  depth,  or  level  flight.  They  can  be  extended  to  a  two 
dimensional  pathkeeping  simulation  by  a  coordinate  transformation.  After  coordinate  rotation 
by  an  angle  P  (positive  in  the  same  direction  as  Q ),  the  resulting  system  is: 


w  =  a,,M>v+ +  0(3  sin(0')  +  +  Fj 


9  =a2,MH'+fl22“9'*‘^23  Sin(0')+^2“^^'^  ^il 


&  =  q 


z'  =  wcos(6')  -  Msin(0') 
x'  =  >vsin(0')  +  Mcos(0') 


(18) 

(19) 

(20) 

(21) 

(22) 


where; 


II 

1 

(23) 

x’=  -z  sin(/J)  +  a:cos(^) 

(24) 

z*  =  zcos(i8)  + Jcsin(/3) 

(25) 

/ 

0,3  =a,3Cos(j3) 

(26) 

/ 

^23  -Clj^COSiP) 

(27) 
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Fj  =  Fj  +a„  cos(0')sin(i3) 
Mj  =  Mj  +023  cos(0')sin(j3) 


(28) 

(29) 


If  the  expected  angular  deviation  from  the  planned  path  is  small,  Equations  (28)  and 
(29)  can  be  simplified  by  assuming  that  cos(6')  is  equal  to  one.  Then  the  rotated  equation  set, 
Equations  (18)  through  (22),  is  identical  in  form  to  Equations  (7)  through  (11). 

Equations  (23)  through  (29)  allow  any  vertical  plane  path  consisting  of  a  series  of 
straight  line  segments  to  be  simulated  one  segment  at  a  time. 

D.  THE  DARPA  SUBOFF 

1.  Background 

For  the  purpose  of  this  work,  it  was  desired  to  have  a  vertical  plane  model  of 
submarine  dynamics  which  would  give  a  similar  response  to  a  modem  fast  attack  nuclear 
submarine  (SSN).  Several  sets  of  unclassified  hydrodynamic  coefficients  were  available,  these 
being  for  the  swimmer  delivery  vehicle  (SDV)  detailed  in  Smith,  Crane,  and  Summey  (1 978) 
and  for  the  DARPA  SUBOFF  model  detailed  in  Roddy  (1990). 

The  SDV  had  a  very  complete  set  of  hydrodynamic  coefficients  which  have  been  used 
in  a  large  number  of  Autonomous  Underwater  Vehicle  (AUV)  research  projects.  Among  these 
is  the  Naval  Postgraduate  School  (NPS)  AUV  sliding  mode  controller,  Hawkinson  (1990). 
Despite  these  advantages,  the  SDV  hydrodynamic  coefficients  were  not  used  because  the  wing 
like  hull  of  the  SDV  bore  little  resemblance  to  an  axisymetric  submarine  hull. 

The  SUBOFF  hydrodynamic  coefficients  detailed  in  Roddy  (1990)  lacked  some  of  the 
cross  coupling  coefficients.  The  documentation  also  lacked  details  on  the  models  metacentric 
height.  Because  the  SUBOFF  represented  a  submarine  hull  form  and  most  of  the  vertical 
plane  coefficients  and  parameters  were  available,  it  was  chosen  as  the  model  for  this  thesis. 

2.  SUBOFF  known  parameters  and  coefficients 

The  SUBOFF  was  developed  to  allow  comparison  between  flow  field  predictions  and 
model  test  data  (Roddy,  1990).  The  available  coefficients  were  based  on  planar  motion 
mechanism  tests  conducted  on  the  model. 
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Because  the  aim  of  the  study  was  to  examine  full  scale  submarine  motions,  the  model 
and  its  hydrodynamic  coefficients  were  scaled  to  a  length  of  300  feet.  After  scaling,  several 
parameters  had  to  be  modified  or  assumed  to  give  control  and  response  comparable  to  a 
modem  fast  attack  submarine.  The  force  coefficients  of  the  stern  planes  was  doubled  to 
provide  a  more  realistic  level  force.  Bow  planes  were  assumed  to  have  one  half  the  force  and 
one  quarter  the  moment  authority  of  the  stern  planes.  Finally,  a  metacentric  height  of  one  foot 
was  assumed,  as  it  provided  a  realistic  point  of  stem  planes  reversal.  The  resulting  parameters 
are  shown  in  Table  1 . 
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Parameter 

SUBOFF  Model 

Scaled  /  Modified  Result 

Length  (Feet) 

14.2917 

300 

Displacement  (tons) 

0.7704 

7,7145 

Maximum  Diameter  (Feet) 

1.667 

35 

Metacentric  Height  (Feet) 

Not  Provided 

1 

XG 

0.00975 

0 

ZG 

Not  Provided 

1 

XB 

-0.006669 

0 

ZB 

Not  Provided 

0 

25, 

-0.005603 

-0.011206 

-0.002409 

-0.004818 

Not  Provided 

-0.005603 

Not  Provided 

0.0012045 

Table  1.  SI 

JBOFF  Assumed  and  modified  parameters 

E.  CONCLUDING  REMARKS 

A  simplified  model  of  submarine  vertical  plane  dynamics  was  derived.  The  coefficients 
for  use  in  this  model  were  obtained  from  the  DARPA  SUBOFF  model,  which  is  a 
representative  axisymetric  submarine  hull  form.  The  simplified  nonlinear  equations  of  motion 
were  incorporated  in  a  SIMULINK®  model  to  allow  easy  integration  with  wave  force  models 
and  different  controllers. 
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Axial  Length  (FI) 


Figure  4.  DARPA  SUBOFF  model,  Roddy  (1990) 
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III.  WAVE  FORCE  MODELING 


A.  INTRODUCTION 

As  a  submarine  operates  near  the  free  surface,  it  encounters  complex  forces  which  may 
cause  unsatisfactory  or  unstable  depth  control.  The  lift  and  moment  from  incident  waves 
increase  in  an  exponential  manner  as  the  surface  is  approached.  To  maintain  a  desired  depth, 
the  ship’s  ballast  is  adjusted  to  counteract  steady  forces.  Control  surfaces  are  used  to  counter 
dynamic  changes.  A  small  depth  excursion  or  change  in  forces  can  overwhelm  the  planes  and 
cause  a  loss  of  depth  control.  The  consequences  range  from  losing  radio  reception  to 
compromising  the  ship’s  mission. 

The  effects  of  incident  waves  on  a  submerged  body  can  be  divided  up  in  several 
categories.  The  largest,  the  first  order  forces  act  at  the  incident  wave  frequency.  These  forces 
move  the  submarine,  but  usually  result  in  oscillations  about  a  mean  state.  Second  order  forces, 
which  are  the  result  of  wave  diffraction  and  wave  interaction,  have  several  different  frequency 
components. 

Wave  diffraction  of  a  single  frequency  wave  results  in  a  steady  force  and  a  varying 
force  at  twice  the  wave  frequency.  The  double  frequency  force  is  typically  neglected,  as  the 
large  inertia  of  the  submarine  effectively  filters  it.  Interactions  of  waves  at  different 
frequencies  also  results  in  forces.  These  consist  of  a  component  acting  at  the  sum  of  the  wave 
frequencies  and  a  component  acting  at  the  difference  of  the  wave  frequencies.  The  sum 
frequency  force  is  typically  neglected,  as  it  is  also  filtered  by  submarine’s  inertia.  The  difference 
frequency  component  results  in  a  slowly  varying  force  on  the  submarine. 

The  slowly  varying  forces  are  the  principle  cause  of  difficult  periscope  depth  control 
(Ni,  Zhang,  and  Dai,  1994).  They  are  compensated  for  using  control  surfaces  and  occasional 
adjustments  of  trim. 

During  the  design  phase,  engineering  decisions  are  made  which  will  determine  the 
ship’s  ability  to  remain  at  periscope  depth.  Of  these,  the  most  critical  are  the  height  of  the  sail 
and  control  surface  sizes.  Every  foot  added  to  the  sail  gives  a  deeper  periscope  depth.  Larger 
planes  improve  the  operator’s  ability  to  compensate  for  changes  in  suction  forces.  However, 
these  improvements  are  not  without  cost.  The  sail  and  other  appendages  are  a  large,  fraction  of 
the  total  drag,  and  can  restrict  the  ship’s  top  speed.  Larger  movable  control  surfaces  can 
adversely  affect  the  high  speed  casualty  recoverability  (Jackson,  1992,  p.  15-9). 


15 


The  goal  of  this  thesis  is  not  to  provide  new  tools  for  the  designer,  but  rather  new 
means  to  enhance  control  for  the  operators  of  current  submarines.  Due  to  this  focus, 
simplified  means  of  modeling  the  wave  forces  for  a  few  specific  cases  will  be  used. 

B.  REVIEW  OF  LINEAR  DEEP  WATER  WAVES 

The  pertinent  features  of  linear  deep  water  waves  will  be  reviewed  to  provide 
background  for  the  following  sections.  The  coordinate  system  used  for  the  examples  is  shown 
in  Figure  5.  For  the  examples  in  this  section,  it  will  be  assumed  that  the  submarine  is  oriented 
with  the  bow  pointing  into  the  page.  Consistent  with  the  global  coordinate  system  from 
Chapter  II,  the  distance  from  the  surface  to  the  submarine  centerline  is  z  .  The  submarine 
diameter  is  D. 


Wove  speed, c 

Hi 

-  L 

d 

Surface  elevation 
shown  at  t  =  0 

Figure  5.  Coordinate  Definition  for  plane  progressive  wave,  adapted  from  Sarpkaya  and 

Isaacson  (1981,  p.  151) 

For  a  wave  of  wavelength  L,  a  wave  number,  k,  can  be  defined. 

^  =  (30) 

L 


Assuming  that  fluid  is  incompressible  and  inviscid  Laplace’s  equation  can  be  applied. 
It  is  thus  desired  to  find  a  solution  to: 


d^(l>  ^ 


(31) 


To  this,  the  boundary  conditions  at  the  free  surface,  and  the  bottom  must  be  applied: 
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—  =  0  ,  at  z  =  d  (no  flow  through  ocean  bottom) 
dz 


(32) 


—  +  +  ^  =  0  ^at  z  =  ri  (zero  velocity  normal  to  ocean  surface) 

dt  dx  dx  dz 


(33) 


9$  1 


+  gri  =  f{t),At  z  =  r) 


(34) 


For  small  amplitude  waves  in  deep  water,  the  following  solution  can  be  obtained 
(adapted  from  Sarpkaya  and  Isaacson,  1981,  p.  159): 

^  TtH  ,  . 

<b  = - e  sin(ajf) 

kT 

=£il  =  l 

k 

L  =  ^ 

2n 

77  =  — cos(<tW) 

^  cos(ct)f) 

^  =  ci)^e~’^  sin(ODf) 

H 

-CxP'  cos{(ot) 

^  sm(cot) 

H 

^  cos{cot) 

H 

^  =  (0^  ~e~^^  sin(£t)f) 


(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 
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where  77  is  the  distance  from  the  surface  to  the  average  level  (z  =  0),  ft)  is  the  angular 
frequency  of  the  incident  wave,  f  is  the  displacement  of  a  particle  in  the  x  direction,  and 
is  the  displacement  of  a  particle  in  the  z  direction. 

A  key  parameter  in  oscillating  flows  is  the  Keulegan-Carpenter  number: 

(^5) 

~D~ 

where  U,^„  is  the  average  velocity  across  the  characteristic  dimension  D. 

By  taking  the  average  of  the  velocity  given  in  Equation  (40) ,  and  substitution  into 
Equation  (45),  the  expression  for  the  Keulegan-Carpenter  can  be  reduced  to  the  following: 


(46) 

D 


Equation  (46)  is  the  Keulegan-Carpenter  number  based  on  the  cross  flow  velocity  of  the 
undisturbed  wave  at  the  same  depth  as  the  centerline  of  the  submarine  hull. 

C.  WAVE  FORCE  REGIMES 

There  are  different  regimes  of  interaction  between  a  submerged  body  and  a  wave  field. 
Broadly,  they  can  be  broken  into  several  areas.  Inertial  interaction,  where  the  body  acts  like  a 
particle  in  the  wave  field.  Wave  diffraction,  where  the  bodies  influence  upon  the  wave  field  is 
accounted  for.  Finally,  there  are  flow  separation  (viscous)  effects.  The  relative  importance  of 
each  of  these  effects  can  be  determined  by  examining  the  relationship  the  body  size  to  the 
wave  parameters.  (Sarpkaya  and  Isaacson,  1981,  pp.  381-386) 
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Figure  6.  Wave  force  regimes  (Sarpkaya  and  Isaacson,  1981,  pg.  385) 


To  estimate  the  significant  effects  for  a  typical  SSN,  a  typical  operating  condition  is 
assumed.  For  a  300  foot  submarine  with  a  35  foot  diameter,  a  typical  periscope  operating 
depth  would  be  about  50  feet  from  the  centerline  of  the  ship  to  the  free  surface.  Using  average 
values  for  sea  states  three  and  four  and  assuming  deep  water  compared  to  the  wavelength,  the 
following  quantities  were  calculated  at  a  depth  of  50  feet: 


Parameter 

Sea  State  3 

Sea  State  4 

Significant  Wave  Height 

3 

6 

Average  Period 

6.623501 

7.154522 

Wave  Length 

224.6467 

262.1114 

Wave  Number 

0.027969 

0.023971 

K 

0.042339 

0.103414 

D/L 

0.1558 

0.133531 

Table  2.  Estimated  Wave  Loading  Parameters 


The  Diameter /Wavelength  (D/L)  ratios  and  the  Keulegan-Carpenter  numbers  of 
Table  2  show  that  for  the  sea  states  of  interest,  wave  diffraction  much  more  significant  than 
viscous  forces.  It  can  be  concluded  that  an  inviscid  analysis  should  give  good  results  for  the 
wave  forces.  However,  this  is  only  rigorous  for  an  unappended  huU,  as  the  control  surfaces 
and  sail  on  an  actual  submarine  will  experience  viscous  effects. 
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D.  SOLUTION  FROM  SLENDER  BODY  THEORY 

Wave  force  solutions  for  several  specific  cases  were  generated  for  the  SUBOFF  by  the 
SSBN  Security  Department  of  the  Johns  Hopkins  University  Applied  Physics  Laboratory.  A 
slender  body  solution  with  some  three  dimensional  corrections  was  used.  The  specific  method 
used  for  the  generation  of  the  first  order  motions  and  second  order  forces  is  detailed  by  O  Dea 
and  Barr  (1976,  pp.  7-25). 

A  seaway  approximation  consisting  of  a  small  number  of  regular  waves  was  used  to 
model  sea  states  three  and  four.  For  each  sea  state,  the  resulting  data  were  separated  into  two 
categories.  The  effects  of  the  first  order  forces  were  given  in  terms  of  body  motions.  The 
effects  of  the  steady  second  order  forces  and  the  difference  interaction  forces  were  provided  in 
pounds  force. 

1,  Seaway  model 

A  random  seaway  can  be  represented  by  the  superposition  of  a  large  number  of  regular 
waves.  The  seaway  was  approximated  by  superimposing  n  regular  waves.  The  frequency  and 
height  of  these  waves  was  determined  using  the  Bretschneider  spectrum.  It  gives  the  spectral 
density  in  terms  of  the  significant  wave  height,  ,  and  the  peak  frequency,  a)„ . 

To  model  sea  state  three,  a  significant  wave  height  of  three  feet  was  used,  with  a  central 
frequency  of  0.836  radians  per  second.  This  results  in  the  following  spectrum: 
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Figure  7  gives  a  statistical  picture  of  the  seaway,  but  is  not  immediately  useful  for  time 
domain  simulation.  One  way  to  obtain  a  time  history  is  to  represent  this  stationary  process  as 
a  the  sum  of  a  series  of  sine  waves; 


n 

tj(0  =  S  A- sin(ft),/ +  a, ) 

/=1 


(48) 


Where  A,  is  the  amplitude  of  the  wave  ,  and  a,  is  its  randomly  chosen  phase  angle. 

If  the  number  of  sine  waves  is  reasonable  large,  and  the  frequencies  and  amplitudes  of 
each  component  are  chosen  to  achieve  the  same  energy  as  the  section  of  spectmm  it 
represents.  Equation  (48)  will  give  a  good  representation  of  the  ocean  surface. 

The  method  chosen  was  to  divide  the  spectra  into  n  segments  of  equal  areas.  This 
results  in  n  sine  waves  all  with  equal  amplitudes.  Integration  of  Equation  (48)  yields: 
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J  S(co)do) 
0 


16 


(49) 


Because  the  spectrum  extends  to  infinity,  it  was  chosen  truncate  the  spectrum  at  a 
point  where  the  area  was  a  fraction  C  of  the  total  area.  The  amount  of  area  to  be  represented 
by  each  sine  wave  is  equal  to  its  mean  square  value.  So  the  amplitude  of  each  sine  wave  is 
equal  to  the  square  root  of  the  area  it  represents  times  the  square  root  of  two. 


A,= 


(50) 


Equation  (48)  can  be  integrated  up  to  some  frequency  ,  which  represents  the 
frequency  at  which  the  spectral  area  is  equal  to  iC  /  n  times  the  total  area. 


(51) 


Solving  Equation  (51)  for  to,  yields; 


0),  =  £0, 


4.  (a 
—In  — 
5  [n 


(52) 


Because  the  spectral  level  is  insignificant  below  a>  equal  to  0.6ct)„ ,  the  frequency  if  the 
first  segment  was  determined  as  follows: 

(0.6a)„  +c5i)  (53) 


CO,  =■ 


The  remainder  of  the  frequencies  were  determined  by  taking  the  midpoint  of  the 
frequencies  at  either  side  of  the  area  segment. 


0):  =  • 


,  for  i  =  2  to  n 
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Figure  8  illustrates  the  method  used,  approximating  the  spectrum  with  sinusoids. 
Nineteen  equal  area  sections  are  divided,  with  the  center  frequency  of  each  segment  marked 
with  a  circle. 


Figure  8.  Spectra  area  division  and  mean  frequencies 

Figure  9  shows  the  ocean  surface  which  results  from  the  use  of  this  method  for  the 
case  of  sea  state  three,  peak  frequency  of  0.862  radians  per  second.  Nineteen  sinusoids  were 
used  to  approximate  the  spectra,  and  the  phase  angles  were  randomly  chosen. 
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2.  First  order  forces 

The  first  order  wave  effects  were  provided  in  the  form  of  submarine  motions.  They 
were  given  as  a  series  of  phasors,  the  real  part  of  the  summation  representing  the  actual 
perturbation  caused  by  the  first  order  wave  forces. 


(55) 

(56) 

Because  the  first  order  motions  were  provided  for  a  specific  depth,  it  was  required  to 
correct  Equations  (55)  and  (56)  for  depth.  The  first  order  motions  roughly  correspond  to  the 
particle  motions  given  by  Equation  (42),  so  an  exponential  decay  was  used  to  derive  the 
following  correction  factor: 
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(57) 


Application  of  Equation  (57)  to  Equations  (55)  and  (56)  results  in: 


(58) 

(59) 

The  displacements  given  by  Equations  (58)  and  (59)  are  not  suitable  for  inclusion  in 
the  submarine  equations  of  motion.  For  this,  an  acceleration  is  required.  Differentiating  twice 
with  respect  to  time  results  in: 


•■(,)  ^  (60) 
6{t)  =  -y'!'  j 

Equations  (60)  and  (61)  were  incorporated  as  force  and  moment  disturbances  in  the 
equations  of  motion  found  in  Chapter  II.  To  test  the  validity  of  this  approach,  an  open  loop 
simulation  was  performed  using  the  accelerations  from  Equations  (60)  and  (61)  for  one  sea 
state  and  heading.  Figure  10  shows  the  results  of  this  simulation,  as  well  as  the  expected  first 
order  motions.  The  upper  curve  shows  the  expected  first  order  motions,  and  the  lower  curve 
shows  the  results  of  integrating  Equations  (7)  through  (1 1)  with  the  accelerations  from 
Equations  (60)  and  (61).  Although  there  was  some  drifting  motion,  the  character  of  motion 
and  the  approximate  amplitude  of  each  cycle  of  motion  very  close.  The  drifting  motion  was  a 
result  lack  of  the  lack  of  open  loop  depth  stability,  which  is  characteristic  of  submarines. 
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Time,  seconds 

Figure  10.  Submarine  response  to  first  order  accelerations,  and  expected  response 

3.  Second  order  forces 

For  a  particular  depth  and  wave  time  history,  the  second  order  forces  were  given  in  the 
following  form; 


(62; 

(63; 

Z(/)  represents  the  force  acting  at  the  body  fixed  coordinate  system  in  the  z  direction 
and  M{t)  is  the  moment  acting  about  the  y  axis.  It  should  be  noted  that  Equations  (63)  and 
(65)  include  the  slowly  varying  forces  (i  j)  and  the  steady  forces  {i  =  j). 


It  can  be  determined  from  analysis  of  the  changes  of  second  order  forces  with  respect 
to  depth  given  by  Crook  (1994,  pp.  61,62)  that  the  steady  forces  with  the  following  exponential 
decay  factor: 

^-2*2  (64) 

The  slowly  varying  order  wave  forces  vary  with  depth  according  to  the  sum  of  the  wave 
numbers: 


(65) 


Application  of  Equations  (64)  and  (65)  to  Equations  (62)  and  (63)  results  in: 


(i((|iUi-ci)y|)/+n,+ay)-(*j+*p(z-z„)) 

(66) 

r"  w  (i((U,-<B  )/+a,+ay)-(*,+*y)(z-z„)) 

Lj=i^iA 

(67) 

The  real  portion  of  Equations  (66)  and  (67)  represents  the  steady  and  slowly  varying 
second  order  wave  forces  acting  on  the  submarine,  with  correction  for  depth. 


4.  Inclusion  of  wave  forces  in  equations  of  motion 

The  first  order  accelerations  and  second  order  forces  had  to  be  combined  to  form  the 
force  and  moment  disturbance  accelerations  for  use  in  the  deeply  submerged  equations  of 
motion  (Equations  (7)  and  (8)). 

'FA>) 

MAt). 


zit) 

Zit) 

— 

0{t) 

+  A/"' 

Mil) 

Incident 

(68) 
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E.  CONCLUDING  REMARKS 

An  elementary  review  of  linear  wave  theory  was  presented.  The  case  of  interest,  a 
submarine  at  periscope  depth,  was  examined  to  determine  the  salient  elements  of  its  interaction 
with  the  incident  waves.  The  parameters  su^ested  that  the  major  features  of  the  incident 
wave  effects  on  the  submarine  could  be  determined  by  using  a  potential  analysis  with  inertial 
and  diffraction  forces  accounted  for. 

The  Bretschneider  spectrum  was  used  to  determine  the  spectral  density  functions  of 
the  sea  states  of  interest.  For  the  purpose  of  time  domain  simulation,  the  spectrum  was 
approximated  by  the  superposition  of  a  number  of  regular  waves  with  randomly  chosen  phase 
angles. 

The  first  order  force  transfer  function  and  second  order  forces  response  amplitude 
operators  were  provided  for  the  SUBOFF  for  a  nominal  speed  and  depth.  Approximate  depth 
scaling  was  introduced  to  allow  use  at  depths  other  than  nominal. 
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IV.  STATE  FEEDBACK  CONTROL  AT  PERISCOPE  DEPTH 


A.  INTRODUCTION 

1.  State  feedback  control 

One  popular  means  of  control  is  to  feed  back  the  system  states  after  the  application  of 
linear  gains.  System  response  of  linear  systems  subjected  to  this  type  of  control  is  predictable, 
and  a  variety  of  tools  are  available  for  control  law  gain  selection. 

The  ship’s  control  party  on  a  submarine  with  conventional  indications  does  not  have 
the  full  state  of  the  ship  to  operate  from.  Although  the  actual  instrumentation  may  vary 
somewhat,  in  general  a  few  analog  indications  are  used  in  conjunction  with  a  digital  depth 
indication.  For  this  reason,  various  levels  of  partial  state  feedback  were  used  to  evaluate  the 
effects  of  missing  indications. 

The  use  of  different  state  feedback  schemes  was  felt  to  be  appropriate  to  model  human 
operators.  The  treatment  of  airplane  pilots  as  a  control  law  “has  come  to  be  recognized  as  a 
quasilinear  element  for  random-appearing  tracking  tasks  related  to  piloting.  At  the  same  time, 
the  pilot  retains  spectacular  nonlinear  gain  changing,  mode  switching,  and  goal  seeking 
precognitive  control  capabilities  as  yet  only  partially  explored.”  (Graham  and  McRuer,  1991,  p. 
1093)  In  this  context,  it  was  assumed  submarine  “pilots”  could  be  treated  in  a  similar  fashion, 
with  feedback  from  each  operating  state  determined  with  linear  gains. 

The  use  of  a  first  order  lag  was  considered  to  model  the  combined  human  and  control 
surface  response  time.  It  was  found  that  reasonable  lag  values  (on  the  order  of  a  half  second) 
had  minimal  effect  on  the  control  response  and  corresponding  submarine  motions.  Because  of 
the  computational  expense,  the  control  response  time  was  neglected. 

State  feedback  control  of  the  linear  system 

X  =  Ax  +  Bu  (ti9) 

where: 


/i  e  91”“'" ,  state  matrix 
fi  €  91"“" ,  control  matrix 
j:  e  91"“' ,  state  vector 
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ue3{^\  control  vector 


can  be  expressed  as: 

u^Kx  (70) 

where: 

a:  6  9?“'"  (71) 

The  system  given  by  Equation  (69)  subject  to  the  control  law  given  by  Equation  (70)  has  the 
following  closed  loop  dynamics  matrix: 

A,  =  iA  +  BK)  (72) 

The  eigenvalues  of  the  closed  loop  dynamics  matrix  will  be  related  to  the  system 
stability  and  responsiveness.  In  general,  the  real  portion  of  the  eigenvalues  must  be  negative 
for  system  stability.  Also  the  more  negative  the  eigenvalues,  the  faster  the  system  response. 

2.  SUBOFF  simulation  parameters 

Wave  force  data  was  available  for  the  SUBOFF  for  four  different  cases.  These  were 
sea  states  three  and  four  with  head  and  beam  directions.  All  were  valid  at  a  speed  of  six  knots 
and  depths  greater  than  fifty  feet. 

At  six  knots,  the  linear  state  representation  used  for  eigenvalue  determination  and 
control  law  design  is: 

0.0179  3.7101  0.0196  OTwl  r-0.0628  -0.10091  (73) 

0.0006  -0.0680  -0.0034  0  q  0.009  -.0027  r^,!  Mjit) 

0  1  0  0  9  ^0  0  0 

1  0  -10.1269  oJLzJ  [  0  0  J  [  0  . 

where: 

(0  ^tritn  ^wave  ftl 
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Mj  (t)  =  +  M^ave(0 


All  simulations  were  performed  using  the  nonlinear  equations  of  submarine  dive  plane 


motion: 

H'  =  a,,«w+a,2«9+ai3Sin(e)+fc,,M^5i+fc,2M^5,  +  F^cos(0)  +  e,,^^+€,29w  (74) 

q  =  a2\uw+a22uq+a23  sm(6)+b2iU^5^  +b22U^5,  +  M j  cos{6)  +  e2\q^  +e22qw  (75) 

e  =  q  (76) 

z  =  w  cos{0)  -  u  sin(6)  (77) 

i:  =  H'sin(0)  +  «cos(6)  (78) 


The  simulations  were  performed  using  a  commanded  depth  of  55  feet  and  using  a  zero 
error  initial  state  vector.  Commanded  pitch  angle,  heave  and  pitch  rate  were  all  zero.  The 
depth  was  chosen  to  provide  a  good  representation  of  actual  submarine  periscope  operating 
depth. 

3.  State  feedback  implementation  with  SIMULINK® 

The  state  feedback  controller  was  implemented  in  the  SIMULINK®  model  shown  in 
Figure  1 1 .  This  block  was  designed  to  use  an  optional  feedforward  signal,  and  also  to  facilitate 
the  use  of  integral  control  (Both  feedforward  and  integral  control  are  discussed  later  in  this 
chapter).  Deflection  limits  were  placed  on  the  control  surfaces.  Control  surface  rate  limits 
were  not  included,  but  could  be  easily  added.  These  limits  are  of  interest  because  of  the 
relationships  between  control  surface  rates,  hydraulic  plant  size  requirements  and  noise  from 
control  surface  operations. 
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A  SIMULJNK®  model  was  developed  to  incorporate  the  submarine  dynamics  of 
Chapter  II,  the  wave  forces  of  Chapter  III,  and  the  state  feedback  control  law.  Also  included 
was  a  logical  means  of  adjusting  the  submarine’s  trim.  This  was  done  by  adding  ballast  in  units 
of  thousands  of  pounds  at  the  center  of  buoyancy,  and  shifting  ballast  from  the  forward  trim 
tank  to  the  after  trim  tank  in  units  of  thousands  of  pounds.  The  details  of  the  trim  model  are 
shown  in  Figure  12,  while  the  overall  model  is  shown  in  Figure  13. 
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Figure  13.  SIMULINK®  state  feedback  control  submarine  model 


4.  Integral  control  on  depth 

To  apply  integral  control,  an  additional  state  is  introduced.  Equations  (74)  through 
(77)  are  augmented  by: 


Z/ 


Z  ^commanded 


(79) 


which  is  used  to  provide  state  feedback.  This  forces  the  steady  state  value  of  ^  to  zero.  In 
general,  this  approach  is  satisfactory  as  long  as  the  control  effort  does  not  become  saturated 
and  the  eigenvalues  of  the  integral  state  are  slower  than  the  state  which  is  being  zeroed. 

5.  Feedforward  of  wave  forces 

Given  the  wave  forces  values,  control  effort  can  be  directly  applied  to  eliminate  the 
average  depth  error.  With  a  constant  disturbance,  a  steady  state  value  of  the  depth  error  can  be 
determined  (Appendix  B).  Using  the  linear  equations  of  motion,  the  steady  state  depth  error 
can  be  written  as  a  linear  combination  of  the  net  force  and  moment  disturbances: 

Z.V.V  ~  ^commanded  ~  ^d  ^ d 
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To  eliminate  the  depth  error,  it  is  required  to  apply  the  control  effort  it  applied: 


K 


24 


(81) 


Equations  (80)  and  (81)  can  be  combined  to  give  a  matrix  gain  relationship  between  the  net 
disturbance  and  the  feedforward: 


14 


24  J 


[C,  Cj] 


F, 


C^K 


2^14 


^1^24  ^2^24 


Fu 


(82) 


The  state  feedback  control  law  with  feedforward  is: 


^bp 


=  Kx  +  Ks 


(83) 


It  has  been  suggested  ( Musker,  Loader,  and  Butcher,  1988),  ( Ni,  Zhang,  and  Dai, 
1994)  that  effective  periscope  depth  control  can  be  achieved  by  feeding  forward  the  average 
second  order  wave  forces.  Because  wave  forces  are  a  dynamic  disturbance  and  the  feedforward 
was  calculated  for  a  steady  disturbance,  a  filter  was  employed  to  cut  out  the  high  frequency 
components  of  the  wave  forces.  The  filter  employed  was  a  first  order  Butterworth  filter  with  a 
cut  off  frequency  .  The  cut  off  frequency  was  initially  chosen  as  one  radian/second.  This 
was  well  below  the  maximum  frequency  wave  force  components  ( around  2.2  radians /second). 
Figure  14  and  Figure  15  show  the  effects  of  the  first  order  Butterworth  filter  on  the  wave 
forces  at  a  depth  of  55  feet  in  sea  state  three.  It  is  apparent  that  with  a  cutoff  frequency  often 
radians  per  second,  the  filtered  forces  and  moments  are  very  close  to  the  unfiltered.  At  the 
lower  cutoff  frequency,  0.1  radians  per  second,  the  filtered  forces  and  moments  are  much 
closer  to  the  average  values. 

To  implement  the  feedforward  control  law,  it  was  assumed  that  the  net  external  force 
and  moment  were  known  quantities.  Equation  (82)  was  implemented  in  the  SIMUUNK 
model  shown  in  Figure  16  while  the  complete  system  model  is  shown  in  Figure  17. 
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Figure  16.  SIMULINK®  model  of  feedforward  calculator 


Figure  17.  SIMULINK®  model  of  system  with  feedforward  term 


6.  Optimization  algorithm  and  parameters 

One  difficulty  of  using  partial  state  feedback  is  that  conventional  pole  placement  or 
linear  quadratic  regulator  algorithms  can  not  be  used  to  determine  the  gains.  The  gains  in 
question  were  selected  randomly,  and  gain  combinations  which  gave  stable  eigenvalues  were 
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simulated.  Because  of  the  clamping  on  the  maximum  planes  angle,  some  gains  which  yielded 
stable  eigenvalues  resulted  in  unstable  ship  control. 

Randomly  selected  gains  certainly  provide  less  than  optimum  depthkeeping.  Because 
of  this,  each  feedback  case  was  optimized  to  provide  the  best  case  for  a  particular  sea  state  and 
commanded  depth  combination.  In  conjunction  with  the  feedback  optimization,  the  trim  was 
optimized. 

The  MATLAB®  function  CONSTR  was  used  to  perform  the  optimizations. 

CONSTR  uses  the  Broyden-Fletcher-Goldfarb-Shanno  variable  metric  method,  and  supports 
constrained  optimization  problems.  To  prevent  the  optimizer  from  selecting  unstable  systems, 
a  constraint  was  placed  on  the  eigenvalues.  The  real  part  of  the  eigenvalues  was  required  to  be 
less  than  a  maximum  value,  usually  -10'^. 

The  objective  of  the  optimizations  was  to  reduce  the  root  mean  square  (RMS)  value  of 
the  depth  error.  For  the  basic  state  feedback  control,  the  average  depth  was  expected  to  differ 
somewhat  from  the  commanded  depth  of  55  feet.  Because  of  this,  the  objective  for  these 
optimizations  was  to  minimize  the  RMS  value  of  the  difference  between  the  depth  and  the 
mean  depth. 

Because  the  optimizations  were  performed  without  regard  for  minimizing  control 
effort  and  or  rates,  large  gains  with  attendant  control  chatter  was  expected.  Although  control 
chatter  is  not  consistent  with  normal  submarine  operations,  it  was  neglected  to  provide  a  clear 
basis  of  comparison  between  the  differing  levels  of  feedback, 

B.  FEEDBACK  OF  DEPTH  AND  PITCH  ANGLE 

1.  Basic  control 

An  elementary  level  of  ship  control  can  be  conducted  with  the  stem  and  bow 
planesman,  each  operating  to  control  one  particular  state.  The  logical  approach  to  this  is  for 
the  stern  planesman  to  control  the  ship’s  angle,  and  the  bow  planesman  to  control  depth.  This 
results  in  the  following  control  law: 


'^commanded 


r^’ 

0 

0 

0 

K,; 

^  ^  commanded 

KJ 

0 

0 

^2.1 

0  _ 

^  ^commanded 

^  ^commanded 


(84) 
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(85) 


|5l<5, 


After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimized  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optimization  statement 
(Vanderplaats,  1984,  p.  9)  is: 


Minimize: 


(86) 


where: 


z  =  depth ,  determined  by  nonlinear  simulation 
'/ 

J  (z)dt 
_  _o _ 

^mean  *  > 

H  =  Ballast  added  to  center  of  buoyang,  thousands  of  pounds 
F  =  Ballast  shiped from  forward  to  aft,  thousands  of pounds 


Subject  to: 


real{eigenvalues{A^))  < 

■^max 


(87) 


Deviation  from  the  mean  value  of  depth,  vice  the  commanded  was  used  because  of  the 
expected  average  depth  error. 

This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  18.  The  results  of  the  four  optimizations  are 
shown  in  Table  3.  For  the  RMS  error  and  maximum  error,  the  optimized  values  are  given, 
along  with  their  percentage  of  the  initial  values. 

In  all  cases,  use  of  the  optimization  resulted  in  reduction  of  the  mean  square  depth 
error  (measured  from  the  average  depth).  Reduction  of  the  maximum  error  was  also  achieved. 
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Sea  State/Direction 

3/head 

3 /beam 

4/head 

4/beam 

Initial  Values 

•••'  •  •  • 

K,4 

0.1465 

0.1465 

0.1465 

0.1465 

K23 

17.51 

17.51 

17.51 

17.51 

H/F  (103  pounds) 

15/0 

15/0 

15/0 

15/0 

Mean  Depth  (feet) 

55.15 

55.20 

55.09 

55.29 

RMS  Error  (feet) 

0.9220 

0.9210 

1.23 

1.30 

Maximum  Error  (feet) 

2.46 

2.47 

3.86 

4.76 

Eigenvalues 

-0.0074  +  0.2096i 

-0.0074  +  0.2096i 

-0.0074  +  0.2096i 

-0.0074  +  0.2096i 

-0.0074  -  0.2096i 

_0.0074  -  0.2096i 

-0.0074  -  0.2096i 

-0.0074  -  0.2096i 

-0.0356  +  0.1144i 

-0.0356  +  0.1144i 

-0.0356  + 0.1 144i 

.  -0.0356  +  0.1144i 

-0.0356  -  0.1144i 

-0.0356- 0.1 144i 

-0.0356 -0.1144i 

-0.0356- 0.1 144i 

Optimized  Values 

' . 

r  ■-  ■  1 

‘i£  •  L!'1 

K,4 

0.567 

0.2293 

0.4708 

0.2016 

Kn 

63.83 

22.5724 

48.6186 

19.58 

H/F  (10^  pounds) 

9.9  /  2.0 

15.8/-4.9 

15.6/-1.5 

4.9/-5 

Mean  Depth  (feet) 

54.7 

55.99 

55.12 

55.44 

RMS  Error  (feet) 

0.4550  (49%) 

0.7549  (82%) 

0.657  (53%) 

1.23  (95%) 

Maximum  Error  (feet) 

1.533  (62%) 

2.03  (82%) 

2.54  (66%) 

4.15  (87%) 

Eigenvalues 

-0.0388  +  0.2392i 

-0.0419  +  0.1464i 

-0.0419  +  0.2178i 

-0.0010  +  0.21 94i 

-0.0388  -  0.2392i 

-0.0419  -  0.1464i 

-0.0419- 0.21 78i 

-0.0010  -  0.21 94i 

-0.0042  +  0.391  li 

-0.0010  +  0.2346i 

-0.0010  +  0.3397i 

-0.0419  +  0.1358i 

-0.0042- 0.391  li 

-0.0010  -  0.2346i 

-0.0010 -0.3397i 

-0.0419  -  0.1358i 

Table  3.  Optimized  pitch  and  depth  control  law  results  and  performance 


delta  rac 


Figure  18.  Simulation  with  depth  and  pitch  angle  control  in  sea  state  three  (head  sea 

direction) 

2.  Disturbance  feedforward 

The  pitch  angle  and  depth  feedback  control  can  be  implemented  with  a  disturbance 
feedforward  to  correct  average  depth  error.  This  results  in  the  following  control  law: 


^  ^commanded 

^bp 

'0 

0 

0 

K,: 

^  ~  ^  commanded 

0 

0 

K2^ 

0  _ 

fi  —  fl 

^  ^commanded 

^  Z  ^commanded  _ 

0 

0 

_M,_ 

where  is  given  by  Equation  (82)  and  the  force  and  moment  disturbances  are  filtered. 
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After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimized  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optimization  statement  is: 


F(K^,,K2^,co,„,H,F)  = 


(90) 


Subject  to: 

real{eigenvalues{A^))  <  (91) 

This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  19.  The  results  of  the  four  optimizations  are 
shown  in  Table  4.  For  the  RMS  error  and  maximum  error,  the  optimized  values  are  given, 
along  with  their  percentage  of  the  initial  values. 
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Sea  State/Direction 

3/head 

3/beam 

4/head 

4/beam 

Initial  Values 

K,4 

0.1465 

0.1465 

0.1465 

0.1465 

Kzi 

17.51 

17.51 

17.51 

17.51 

1 

1 

1 

1 

H/F  (10^  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.07 

56.7 

55.826 

61.33 

RMS  Error  (feet) 

0.408 

2.24 

1.71 

6.93 

Maximum  Error  (feet) 

1.104 

5.27 

3.71 

13.46 

Eigenvalues 

-0.0074  +  0.2096i 

-0.0074  +  0.20961 

-0.0074  +  0.20961 

-0.0074  +  0.20961 

-0.0074  -  0.20961 

-0.0074  -  0.20961 

-0.0074  -  0.20961 

-0.0074  -  0.20961 

-0.0356  +  0.11441 

-0.0356  +  0.11441 

-0.0356  +  0.11441 

-0.0356  +  0.11441 

-0.0356-0.11441 

-0.0356-0.11441 

-0.0356  -  0.11441 

-0.0356-0.11441 

Optimized  Values 

K,4 

1.116 

3.5763 

7.40 

3.396 

Kzi 

151.00 

454.7  1 

1,073.6 

979.02 

(rad/sec) 

0.743 

3.30  1 

6.83 

6.43 

H/F  (1  a’ pounds) 

19.5/3.5 

22.1/1.3 

26.5/3.25 

8.4/-4.0 

Mean  Depth  (feet) 

54.996 

55.14 

55.04 

55.21 

RMS  Error  (feet) 

0.102  (25%) 

0.556  (25%) 

0.810  (47%) 

0.883  (13%) 

Maximum  Error  (feet) 

0.322  (29%) 

2.24  (43%) 

0.560  (15%) 

3.36  (25%) 

Eigenvalues 

-0.0337  +  0.33451 

-0.0354  +  0.60551 

-0.0322  +  0.86041 

-0.0250  +  0.56751 

-0.0337  -  0.33451 

-0.0354  -  0.60551 

-0.0322  -  0.86041 

-0.0250  -  0.56751 

-0.0092  +  0.60881 

-0.0076  +  1.04901 

-0.0107  +  1.62961 

-0.0179  +  1.60191 

-0.0092-0.60881 

-0.0076  -  1.04901 

-0.0107  - 1.62961 

-0.0179  - 1.60191 

Table  4.  Optimized  pitch  and  depth  control  law  with  disturbance  feedforward  results  and 

performance 
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Figure  19.  Simulation  with  depth  and  pitch  angle  control  with  disturbance  feedforward,  sea 

state  three  (head  seas) 

3.  Integral  control 

The  feedback  of  depth  and  pitch  angle  can  be  augmented  with  integral  control  on 
depth  to  remove  the  average  depth  error.  Since  the  bow  planes  are  principally  used  for  the 
control  of  depth,  the  integral  control  was  applied  to  the  bow  planes  only.  This  results  in  the 
following  control  law: 
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^  ^commanded 


0 

0 

0 

^,4 

1 - 

0 

0 

to 

0 

0  J 

^  ^  commanded 

^ 

^  ^commanded 
^  ^commanded 
2/ 


|5|<5, 


(92) 


(93) 


After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimized  to 
minimize  the  deviation  from  the  commanded  depth.  The  formal  optimization  statement  is: 


Minimize: 


F( A"i4  ,  Fi5  ,  ^23  ,H,F)  - 


j 

J  Z commanded  ) 


(94) 


Subject  to: 


real{eigenvalues{A^))  <  F, 


(95) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  20.  The  results  of  the  four  optimizations  are 
shown  in  Table  5. 
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Maximum  Error  (feet) 


Eigenvalues 


-0.0149  +  0.8824i  I  -0.0134  +  0.5221i  I  -0.0001  +  0.261 5i  I  -0.0163  +  0.4296i 


-0.0149  -  0.8824i 
-0.0274  +  0.3887i 
-0.0274  -  0.3887i 


-0.0134 -0.5221i 


-0.0001-  0.261 5i 


-0.0163  -  0.4296i 


-0.0286  +  0.2475i  -0.0421  +  0.1675i  -0.0258  +  0.1 707i 


-0.0286  -  0.2475i 
-0.0020 


-0.0421  -01 675i  -0.0258  -  0.1 707i 


-0.0015 


-0.0015 


Table  5.  Optimized  pitch  and  depth  integral  control  law  results  and  performance 
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delta  rad 


Figure  20.  Simulation  with  depth,  pitch  angle,  and  integral  control,  sea  state  three  (head 

seas) 


C.  FULL  STATE  FEEDBACK  WITH  PARTIAL  DISTRIBUTION 
1.  Basic  control 

The  poor  depth  control  of  the  previous  section  can  be  improved  be  adding  to  the 
number  of  states  fed  back.  In  keeping  with  previous  logic,  the  bow  planes  wiU  be  controUed 
by  the  depth  and  heave,  while  the  stern  planes  will  be  controlled  by  pitch  and  pitch  rate.  This 
results  in  the  following  control  law: 
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A/ 

KJ 

0  0 
^22  ^23 


Ku 

0 


^  ^commanded 

^  Q  commanded 
n  _ 

^  '^commanded 
^  ~  ^commanded  _ 


15|<5, 


(96) 


(97) 


After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimized  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optimization  statement  is: 


Minimize: 


F(K,,,K„,K22,K2^.H,F) 


'/ 


(98) 


Subject  to: 

realieigenvalues{A^))^E^„  (99) 

This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  21.  The  results  of  the  four  optimizations  are 
shown  in  Table  6. 
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Sea 

3/beam 

4/head 

4/beam 

State/Direction 

Initial  Values 

KJ 

12.5543 

0 

12.5543 

0 

12.5543 

0 

12.5543 

0 

0 

22.5268 

0 

22.5268 

0 

22.5268 

0 

22.5268 

0 

24.9900 

0 

24.9900 

0 

24.9900 

0 

24.9900 

2.5490 

0 

2.5490 

0 

2.5490 

0 

2.5490 

0 

H/F  (103  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.06 

55.33 

55.21 

55.73 

RMS  Error  (feet) 

0.216 

0.425 

0.412 

1.06 

Maximum  Error  (feet) 

0,796 

1.38 

1.26 

3.38 

Eigenvalues 

-0.6743 

-0.6743 

-0.6743 

-0.6743 

-0.0413  +  0.3587i 

-0.0413  +  0.35871 

-0.0413  +  0.35871 

-0.0413  +  0.35871 

-0.0413 

-  0.3587i 

-0.0413 

-  0.3587i 

-0.0413 

-  0.3587i 

-0.0413 

-  0.3587i 

-0.1789 

-0.1789 

-0.1789 

-0.1789 

Optimized  Values 

KF 

0 

7.035 

0 

13.265 

0 

16.58 

0 

0 

1366.3 

0 

91.224 

0 

91.49 

0 

96.44 

0 

1163.3 

0 

20.87 

0 

51.15 

0 

82.43 

89.26 

0 

3.813 

0 

2.1048 

0 

15.33 

0 

H/F  (103  pounds) 

14.6/-1.4 

12.3/-0.4 

20.0/0.1 

16.8/0.0 

Mean  Depth  (feet) 

55.02 

55.01 

55.21 

55.09 

RMS  Error  (feet) 

0.1969  (91%) 

0.350  (82%) 

0.358  (87%) 

0.821  (77%) 

Maximum  Error  (feet) 

1.17  (147%) 

0.99  (72%) 

1.13  (90%) 

3.46  (102%) 

Eigenvalues 

-3.3197 

-0.1266  +  0.47091 

-0.7286 

-0.2213  +  0.88611 

-0.0693  +  1.25221 

-0.1266 

-  0.4709i 

-0.1458  +  0.47001 

-0.2213 

-  0.8861i 

-0.0693 

-  1.2522i 

-0.2619  +  0.01791 

-0.1458 

-  0.4700i 

-0.3797  +  0.48191 

-0.1879 

-0.2619 

-  0.01 79i 

-0.1514 

-0.3797 

-0.48191 

Table  6.  Full  state  feedback  (partial  distribution)  control  law  optimization  results  and 

performance 
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delta  rad  2  feet 


Figure  21.  Simulation  with  fuU  state  partial  distribution  feedback  control,  sea  state  three 

(head  seas) 


2.  Disturbance  feedforward 

As  before,  the  basic  control  law  can  be  modified  to  include  a  feedforward  term  to 
correct  the  average  depth  error. 


^  ^commanded 

■^1. 

0 

0 

^  ~  ^  commanded 

C,K,,- 

'Fu' 

0 

K22 

^23 

0 

Q  __  f) 

^  ^commanded 

Z~‘  Z  commanded  _ 

0 

0 

_Mj_ 

(101) 
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After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimized  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optimization  statement  is: 


Minimize: 


3 

J  “  ^cummanded  ) 


(102) 


Subject  to: 


real{eigenvalues(Ac))  ^  E, 


(103) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  22.  The  results  of  the  four  optimizations  are 
shown  in  Table  7. 
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Sea 

State/Direction 


Initial  Values 


KT 

12.5543 

0 

0 

22.5268 

0 

24.9900 

2.5490 

0 

CO^,,  (rad/sec) 

1 

H/F  (103  pounds) 

20/0 

Mean  Depth  (feet) 

55.16 

RMS  Error  (feet) 

0.371 

Maximum  Error  (feet) 

1.264 

Eigenvalues 

-0.6743 

-0.0413  +  0.3587i 

-0.0413  -  0.3587i 

-0.1789 

Optimized  Values 

KT 

25.81 

0 

0 

175.77 

0 

26.1 

3.5847 

0 

(rad/sec) 

0.481 

H/F  (10^  pounds) 

7.1 /-3.3 

Mean  Depth  (feet) 

55.01 

RMS  Error  (feet) 

0.0785  (21%) 

Maximum  Error  (feet) 

0.246  (19%) 

Eigenvalues 

-1.1709 

-0.5214 

-0.0948 

-0.3993 

12.5543 


12.5543 


12.5543 


2.5490 


0  22.5268 

0  24.9900 

)0  0 


2.5490 


0  22.5268 

0  24.9900 

)0  0 


2.5490 


0  22.5268 

0  24.9900 

)0  0 


20/0 


55.35 


0.533 


1.83 


-0.6743 


20/0 


55.23 


0.551 


1.72 


-0.6743 


1 


20/0 


55.82 


1.40 


3.50 


-0.6743 


-0.0413  +  0.35871  -0.0413  +  0.35871  -0.0413  +  0.35871 


-0.0413  -  0.35871 
-0.1789 


15.73  0 

0  230.65 

0  535.98 


-0.0413  -  0.35871 
-0.1789 


12.74  0 

0  18.09 


-0.0413  -  0.35871 
-0.1789 


16.87  I  0 


0  18.12 


3.80 

0.998 

0.999 

19.7/-2.5 

20.0/0.0 

19.7/0,0 

55.23 

55.18 

55.63 

0.3624  (68%) 

0.5171  (94%) 

1.25  (89%) 

1.07  (58%) 

2.03  (118%) 

3.26  (93%) 

-0.1692  +  1.3601i 

-0.6730 

-0.9159 

-0.1692-1.36011 

-0.0429  +  0.3320i 

-0.0386  +  0.3173i 

-0.6826  +  0.4529i 

-0.0429  -  0.3320i 

-0.0386 -0.3173i 

-0.6826  -  0.4529i 

-0.1767 

-0.1769 

Table  7.  Full  state  feedback  (partial  distribution)  with  disturbance  feedforward  control  law 

optimization  results  and  performance 


51 


delta  rad 


Figure  22.  Simulation  with  full  state  partial  distribution  control  and  disturbance 
feedforward,  sea  state  three  (head  seas) 


3.  Integral  Control 

This  full  state  feedback  with  partial  distribution  was  augmented  with  integral  control 
on  depth  to  remove  the  average  depth  error.  As  before,  the  integral  control  was  done  using 
the  bow  planes  only.  This  results  in  the  following  control  law: 


0  0 

0  K22  ^23  ® 


(104) 


|5|<5, 


(105) 
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After  a  stable  set  of  random  gains  was  determined,  the  controller  was  optimi2ed  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optimization  statement  is  as 
follows: 


Minimize: 


(106) 


Subject  to: 


real{eigenvalues{A^))  < 


(107) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  23.  The  results  of  the  four  optimizations  are 
shown  in  Table  8. 
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Sea  State /Direction 


Initial  Values 


KJ 


H/F  (105  pounds) 


Mean  Depth  (feet) 


RMS  Error  (feet) 


Maximum  Error  (feet) 


3/head 


3/beam 


4/head 


4/beam 


Optimized  Values 


12.5543  0 

0  22.5268 

0  24.9900 

2.5490  0 

0.0010  0 

20/0 


55.05 


0.2106 


0.866 


-0.6744 

-0.0412  +  0.3587i 
-0.0412  -  0.3587i 
-0.1785 
-0.0004 


523.12 

0  1: 

0  i: 

89.26 
.04016 


12.5543  0 

0  22.5268 

0  24.9900 

2.5490  0 

0.0010  0 

20/0 


55.31 


0.533 


1.723 


-0.6744 

-0.0412  +  0.3587i 
-0.0412-0.35871 
-0.1785 
-0.0004 


12.5543 

0 

12.5543 

0 

0 

22.526 

0 

22.526 

0 

8 

0 

8 

2.5490 

24.990 

2.5490 

24.99 

0.0010 

0 

0.0010 

0 

20/0 


55.21 


0.4868 


1.82 


-0.6744 

-0.0412  +  0.3587i 
-0.0412  -  0.3587i 
-0.1785 
-0.0004 


385.2 

0 

0 

254.4 

0.0077 


0 

874.3 

1400.5 

0 

0 


14.074 

0 

0 

2.146 

0.0003 


20/0 


55.74 


1.205 


3.43 


-0.6744 
-0.0412  + 
0.3587i 

-0.0412  -  0.3587i 
-0.1785 
-0.0004 


79.08  I  0 
0  353.83 

0  54.45 

52.88  0 

0.0032 


H/F  (10^  pounds) 

14.6/-1.4 

26.7/4.6 

21.2/-0.2 

10.5/-2.6 

Mean  Depth  (feet) 

55.02 

55.05 

55.25 

55.18 

RMS  Error  (feet) 

0.059  (28%) 

0.3017  (57%) 

0.4274  (88%) 

0.909  (75%) 

Maximum  Error  (feet) 

0.248  (29%) 

0.862  (50%) 

1.16  (64%) 

3.59  (105%) 

Eigenvalues 

-30.7579 

-22.3751 

-0.7495 

-3.4631 

-4.6840 

-1.8146  +  1.6555i 

-0.1526  +  0.3297i 

-1.4563 

-1.0562 

-1.8146-  1.6555i 

-0.1526  -  0.3297i 

-0.9504 

-0.1695 

-0.6572 

-0.1168 

-0.1490 

-0.0005 

0.00003 

-0.0001 

-0.0001 

Table  8.  Full  state  feedback  (partial  distribution)  integral  control  law  optimization  results  and 

performance 
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delta  rad  z  feet 


D.  FULL  STATE  FEEDBACK 
1.  Basic  control 

The  best  control  possible  using  state  feedback  should  result  from  the  use  of  all  four 
states  by  each  control.  This  results  in  the  foUowing  control  law: 
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(108) 


^1. 

K,2 

K,2 

^  ^commanded 

Q  ~  ^  commanded 

^2, 

K22 

^23 

to 

1 _ 

ft  —  0 

^  '^commanded 

_  ^commanded  _ 

(109) 


After  a  stable  set  of  gains  was  determined  using  a  linear  quadratic  regulator  algorithm, 
the  controller  was  optimi2ed  to  minimize  the  deviation  from  the  average  depth.  The  formal 
optimization  statement  was: 

Minimize: 


(110) 


Subject  to: 


real{eigenvalues{A^))  < 


(111) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimum  response  is  shown  in  Figure  24.  The  results  of  the  four  optimizations  are 
shown  in  Table  9. 
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Sea  State/Direction 


Initial  Values 


K'T 


3/head 


4/head 


4/beam 


5.1622 

6.847 

5.1622 

6.847 

-S' '*^■45'.' 

5.1622 

6.847 

12132 

-168.26 

121.32 

-168.26 

121.32 

-168.26 

27.744 

-65.795 

27.744 

-65.795 

27.744 

-65.795 

-0.0789 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.9740 

H/F  (10-^  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.09 

55.37 

55.21 

56.16 

RMS  Error  (feet) 

0.0914 

0.3605 

0.355 

1.60 

Maximum  Error  (feet) 

0.262 

1.34 

1.09 

4.55 

Eigenvalues 

-0.8859 

-0.8859 

-0.8859 

-0.8859 

-0.2854 

-0.2854 

-0.2854 

-0.2854 

-0.1630  +  0.2247i 

-0.1630  +  0.22471 

-0.1630  +  0.22471 

-0.1630  +  0.22471 

-0.1630  -  0.22471 

-0.1630  -  0.2247i 

-0.1630 

-  0.2247i 

-0.1630  -  0.2247i 

Optimized  Values 

KJ 

10.300  9.638 

4.591 

11.31 

4.591 

11.306 

3.503 

11.5734 

45.790  170.45 

-283.9 

78.91 

-283.9 

78.91 

-374.02 

178.92 

-135.16  39.390 

125.7 

-1.333 

-125.7 

-1.333 

-66.65 

40.760 

3.6389  0.0308 

1.457 

0.1568 

1.457 

0.1568 

0.446 

0.0221 

H/F  (10^  pounds) 

20.1 /-0.9 

18.0/-0.3 

18.0/-0.33 

19.9/0.0 

Mean  Depth  (feet) 

55.03 

55.13 

55.13 

56.54 

RMS  Error  (feet) 

0.037  (40%) 

0.2638  (73%) 

0.2683  (76%) 

1.24  (78%) 

Maximum  Error  (feet) 

0.119  (45%) 

0.961  (72%) 

0.961  (88%) 

4.06  (89%) 

Eigenvalues 

-1.7613 

-1.0446 

-1.0446 

-1.0313  +  0.42181 

-0.2510 

-0.3493  +  0.2881i 

-0.3493  +  0.2881i 

-1.0313  -  0.42181 

-0.0617  +  0.52641 

-0.3493  -  0.2881i 

-0.3493 

-  0.2881i 

-0.0936  +  0.0741i 

-0.0617-0.52641 

-0.2053 

-0.2053 

-0.0936  -  0.0741i 
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After  a  stable  set  of  gains  was  determined  using  a  linear  quadratic  regulator  algorithm, 
the  controller  was  optimized  to  minimize  the  deviation  from  the  average  depth.  The  formal 
optimization  statement  was: 

Minimize: 


FCATii  ,  K^2  ’  ^13  »  ^14  ’  ^21  ’  ^22  »  ^23  »  ^24  y  ^ ,  F)  — 


U 


(114) 


Subject  to: 


real{eigenvalues{A^)')  < 


(115) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimum  response  is  shown  in  Figure  25.  The  results  of  the  four  optimi2ations  are 
shown  in  Table  10. 
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Sea  State/Direction 

3/head 

3/beam 

4/head 

4/beam 

Initial  Values 

IMS 

ICT 

6.847 

5.1622 

6.847 

5.1622 

6.847 

5.1622 

6.847 

5.1622 

-168.26 

121.32 

-168.26 

121.32 

-168.26 

121.32 

-168.26 

121.32 

-65.795 

27.744 

-65.795 

27.744 

-65.795 

27.744 

-65.795 

27.744 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.9740 

-0.0789 

(rad/sec) 

1 

1 

1 

1 

H/F  (103  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.06 

55.937 

55.21 

56.70 

RMS  Error  (feet) 

0.1428 

1.343 

0.430 

2.393 

Maximum  Error  (feet) 

0.4897 

3.46 

1.11 

5.49 

Eigenvalues 

-0.8859 

-0.8859 

-0.8859 

-0.8859 

-0.2854 

-0.2854 

-0.2854 

-0.2854 

-0.1630  4-  0.2247i 

-0.1630  +  0.2247i 

-0.1630  +  0.2247i 

-0.1630  +  0.2247i 

-0.1630 

-  0.2247i 

-0.1630 

-  0.2247i 

-0.1630 

-  0.2247i 

-0.1630 

-  0.2247i 

Optimized  Values 

KTT 

-26.2 

96.5 

0.3964 

11.08 

20.60 

-2.681 

5.00 

-3.94 

-1633 

914.4 

-879.14 

329.16 

-194.6 

181.5 

-452.53 

609.1 

-368.8 

91.9 

-297.0 

-91.096 

27.66 

48.6 

25.47 

449.5 

5 

-1.2 

3.52 

0.215 

-0.276 

-0.238 

6.66 

0.318 

CO^„  (rad/sec) 

1.4046 

1.033 

0.983 

1.739 

H/F  (103  pounds) 

13.5/-0.8 

19.16/-0,2134 

19.0/-0.1 

18.8/-0.1 

Mean  Depth  (feet) 

55.0013 

55.18 

55.18 

55.22 

RMS  Error  (feet) 

0.0928  (65%) 

0.4121  (31%) 

0.400  (36%) 

0.792  (33%) 

Maximum  Error  (feet) 

0.2852  (58%) 

1.62  (47%) 

1.117  (101%) 

2.23  (41%) 

Eigenvalues 

-7.0391 

-0.8095  +  0.5268i 

-1.1927 

-0.0898  +  0.92971 

-4.6280 

-0.8095 

-  0.5268i 

-0.2499  +  0.3312i 

-0.0898 

-  0.9297i 

-0.1322  +  0.1329i 

-1.1529 

-0.2499 

-  0.3312i 

-1.1834 

-0.1322 

-  0.1329i 

-0.0315 

-0.0618 

-0.6517 

Table  10.  Full  state  feedback  control  law  with  disturbance  feedforward  optimization  results 

and  performance 
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Figure  25.  Full  state  feedback  control  with  disturbance  feedforward  optimized  control 

simulation,  sea  state  three  (head  seas) 


3.  Integral  control 

This  full  state  feedback  with  partial  distribution  was  augmented  with  integral  control 
on  depth  to  remove  the  average  depth  error.  Since  the  bow  planes  are  principally  used  for  the 
control  of  depth,  the  integral  control  was  done  using  the  bow  planes  only.  This  results  in  the 
following  control  law: 


^  ^commanded 
^  ~~  ^commanded 

ft  —  d 

^  ^commanded 
^  ^commanded 


(116) 


|5l<5. 


(117) 
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After  a  stable  set  of  random  gains  was  determined,  the  controllei  was  optimized  to 
minimize  the  deviation  from  the  average  depth.  The  formal  optirntzation  statement  is: 


Minimize: 


) 

J  (2  ^commanded  ^ 


(118) 


Subject  to: 


real{eigenvalues{Ac))  ^  ^max 


(119) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  opdmum  response  is  shown  in  Figure  26.  The  results  of  the  four  opdmizadons  are 

shown  in  Table  1 1 . 
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Sea  State/Direction 

3/head 

3/beam 

4/head 

4/beam 

Initial  Values 

liW^ 

6.847 

5.1622 

6.847 

5.1622 

6.847 

5.1622 

6.847 

5.1622 

-168.26 

121.32 

-168.26 

121.32 

-168.26 

121.32 

-168.26 

121.32 

-65.795 

27.744 

-65.795 

27.744 

-65.795 

27.744 

-65.795 

27.744 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.9740 

-0.0789 

0.01 

-0.01 

0.01 

-0.01 

0.01 

-0.01 

H/F  (lO-"^  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.01 

55.02 

54.78 

54.68 

RMS  Error  (feet) 

0.101 

0.415 

0.573 

2.483 

Maximum  Error  (feet) 

0.3065 

1.545 

1.79 

6.927 

Eigenvalues 

-0.8854 

-0.8854 

-0.8854 

-0.8854 

-0.2693 

-0.2693 

-0.2693 

-0.2693 

-0.1652  +  0.2153i 

-0.1652  + 0.21 53i 

-0.1652  +  0.2153i 

-0.1652  +  0.2153i 

-0.1652 

-0.2153i 

-0.1652 

-  0.21 53i 

-0.1652 

-  0.21 53i 

-0.1652 

-  0.2153i 

-0.122 

-0.122 

-0.122 

-0.122 

Optimized  Values 

KJ 

240.17 

24.2736 

1.3875 

7.5561 

13.059 

4.681 

-1.679 

6.870 

-137.3 

256.28 

-158.91 

153.309 

-146.54 

155.09 

-234.17 

257.016 

495.76 

84.99 

-52.850 

39.280 

-36.04 

42.31 

-94.42 

48.260 

36.65 

0.1753 

0.4948 

-0.0550 

0.9425 

-0.0484 

1.140 

-0.0327 

0.162 

-0.0128 

0.0109 

0.0069 

0.013 

-0.0037 

0.0111 

-0.0105 

H/F  (10^  pounds) 

20.7/1.9 

19.1/0 

20/0 

19.9/0.0 

Mean  Depth  (feet) 

55.00 

55.00 

54.99 

54.92 

RMS  Error  (feet) 

0.0414  (14%) 

0.372  (90%) 

0.536  (93%) 

1.96  (79%) 

Maximum  Error  (feet) 

0.175  (57%) 

1.01  (65%) 

1.57  (88%) 

6.88  (99%) 

Eigenvalues 

-17.3244 

-0.8739 

-1.2707 

-1.0690 

-0.7017 

-0.1541  +0.1441i 

-0.2086  +  0.1628i 

-0.1214  +  0.3327i 

-0.1942  +  0.4296i 

-0.1541 

-  0.1441i 

-0.2086 

-  0.1628i 

-0.1214 

-  0.3327i 

-0.1942 

-  0.4296i 

-0.2567 

-0.2056 

-0.2428 

-0.0076 

-0.0375 

-0.0208 

-0.0012 

Table  11.  Full  state  feedback  integral  control  law  optimization  results  and  performance 
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(head  seas) 


E.  CONCLUDING  REMARKS 

For  each  case  of  feedback  control,  the  degree  of  control  achieved  generaUy  improved 
by  the  additional  state  feedbacks.  Full  distribution  of  each  state  to  both  controls  fbrther 
reduced  the  error.  Table  12  provides  a  summary  of  the  optimizations  performed,  and  the  RMS 

error  of  each  one. 

Changes  in  the  optimized  trim  and  control  law  in  all  cases  varied  substantially  with 
changes  in  sea  state  ot  direction.  This  is  consistent  with  operational  experience. 

Each  controller  was  optimized  with  only  the  goal  of  minimizing  the  mean  square  of 
the  depth  error.  This  resulted  in  large  gains  and  excessive  control  effort.  In  addition,  large 
rates  of  control  were  experienced.  This  would  be  detrimental  for  actual  submarine  operauons, 
as  there  are  rate  limns  associated  with  the  connol  surfaces.  These  limits  come  from  the  sizing 
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of  the  hydraulic  plants  which  drive  the  planes,  and  operation  concerns  related  to  plane  induced 
noise. 

Some  improvements  in  depthkeeping  were  achieved  by  the  feedforward  of  the 
disturbance  forces.  This  is  in  spite  of  the  feedforward  being  based  on  a  steady  state  response 
to  a  constant  disturbance. 


Sea  State/Direction 

3/Head 

3/ Beam 

4/Head 

4/Beam 

Control  Scheme 

Depth  and  Pitch  Angle 

0.4550 

0.7549 

0.657 

1.23 

Depth  and  Pitch  angle  with  feedforward 

0.102 

0.556 

0.810 

0.883 

Depth  and  Pitch  angle  with  integral 

0.455 

0.3811 

0.865 

1.05 

Full  State  with  partial  distribution 

0.1969 

0.350 

0.358 

0.821 

Full  State  with  partial  distribution  and  feedforward 

0.0785 

0.3624 

0.5171 

1.25 

Full  State  with  partial  distribution  and  integral 

'  0.059 

0.3017 

0.4274 

0.909 

Full  State 

0.037 

0.2638 

0.2683 

1.24 

Full  State  with  feedforward 

0.0928 

0.4121 

0.400 

0.792 

Full  State  integral 

0.0414 

0.372 

0.536 

1.96 

Table  12.  Optimized  RMS  error  (feet)  of  state  feedback  control  schemes 
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V. 


SLIDING  MODE  CONTROL 


A.  INTRODUCTION 

1.  Overview  of  MIMO  sliding  mode  control 

The  controller  design  starts  with  a  standard  linear  state  space  representation; 

x=  Ax  +  Bu  0  20) 


where: 

A  £  ,  state  matrix 

5  e  91"“" ,  control  matrix 
X  £  91""'' ,  state  vector 
M  e  91'“' ,  control  vector 


The  sliding  mode  control  law,  u,  is  composed  of  two  main  parts: 

u  =  u  +  u 


(121) 


The  first  part,  ii ,  is  a  linear  feedback  based  on  the  linear  representation  given  by  Equation 
(120).  The  second  part,  u  ,  are  nonlinear  feedbacks  with  their  signs  switching  depending  on 
the  relationship  of  the  system  states  to  the  sliding  surfaces.  The  sliding  surfaces  are  hyper 
planes  in  the  state  space,  with  one  for  each  control.  They  are  defined  by: 

a(x)  =  s'^x  =  0  (122) 


Where: 


S  £  91"-'''" 


To  determine  the  nonlinear  feedback  functions,  the  concept  of  Liapunov  stability  is 
used.  The  Liapunov  function  is  taken  as: 


V'W=-j((Tf  +(yl+...+al) 


(123) 
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Asymptotic  system  s 
that  is: 


tability  is  guaranteed  provided  that  V(x)  is  a  positive  definite  function. 


V(X)  =  CT|C7,  +  d20'2+-  +^n^n  <  ® 


(124) 


Equation  (124)  is  satisfied  if: 


d.  cT,  <  0  for  i  —1  to  n 


(125) 


Equation  (125)  can  be  rewritten  as: 


dj  =-riiSigni(7j) 


(126) 


where  is  a  posiove  gain  parameter  for  the  sKdtng  surface.  Equadon  (126)  can  be  rewritten 
after  subsrituting  the  rime  derivative  of  Equation  (122)  and  Equation  (120)  as: 


S'^iAx  +  Bu)  =  - 


7],5Ign((T,) 

r]2sign{02) 


r\„sign{a„) 


Solving  Equation  (127)  for  u  yields: 


T]jsign(a,) 

Ti2sign{a2) 


u  =  -(S''Br'S^Ax-(S^B)-' 


ri„sign{o„) 


(128) 


Equation  (128)  can  be  rewritten  in  matrix  form  as: 

«=  Kx+  Kjisign{o) 


(129) 


where: 
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K  = -(S^  By' s'^  A 
K,=-{S'^By' 

Equation  (129)  is  identical  in  form  to  Equation  (121),  with  a  linear  state  feedback  and  a 
nonlinear  switching  term. 

For  the  decomposition  in  Equation  (128),  it  is  required  that  the  closed  loop  stability 
matrix,  (A-BK),  have  n  zero  eigenvalues.  The  sliding  surfaces  are  the  left  eigenvectors  resulting 
from  the  zero  eigenvalues. 

2.  Utkin’s  method  for  MIMO  sliding  mode  control  law  design 
Determination  of  the  sliding  surfaces  can  be  difficult,  especially  with  a  MIMO  system. 
One  technique  for  this  is  proposed  by  Utkin  (1977).  For  this  technique  to  be  applied,  the  B 
matrix  of  the  state  space  system  (Equation  (120))  must  be  of  the  form: 

(130) 

where: 

5,  € 

For  the  MEMO  cases  of  vertical  plane  depth  control  used  in  this  thesis,  this  was  the  case.  For 
the  stern  planes  only  control  examples,  a  QR  factorization  was  applied,  to  transform  the  state 
space  system  into  this  form. 

Given  that  the  B  matrix  is  of  the  form  defined  in  Equation  (130),  Equation  (120)  can 
be  decomposed  into  the  following: 

x,  =  A,|jr|  +  ;4|2a:2  +  S|M  (131) 

^2  —  ^21-^1  ^22^2  (132) 
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where: 


X,  e 

X2  e  gt'™-"’-'' 

A,,  €91"” 

A, 2  e  g?""*"'-"' 

^2,  e  9?^'”'"'” 

^22  6 

The  sliding  surfaces  become: 

a  =  s[  X,  +  5^X2  =  0 

where: 

si  e9t'"" 
si  e9i'“<'"""> 


Because  the  sliding  surfaces  are  the  left  eigenvectors  of  the  «  zero  eigenvalues,  Sl  can 
be  set  to  the  identity  matrix  without  loss  of  generality.  Substitution  of  Equation  (133)  into 
Equation  (148)  leads  to: 

(j  =  X,  +  5^X2  =  -rjsign{X]  +5^X2)  0  ^4) 


T},«g/J(CT|) 

r]2sign(a2) 


u  =  -Bj  '[(4||  +5J 42i)X|  +  (4,2  +  si A22)X2\  -  ' 


ri„sign(a„) 


(135) 
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When  the  system  is  on  the  sliding  surfaces.  Equation  (1 33)  can  be  used  to  solve  for  jt, 
in  terms  of  ,  resulting  in: 


=  -Sixj 


(136) 


Substitution  of  this  result  into  Equation  (132)  results  in: 

•^2  “  (^^22  ~  ^21  *5^  )X2 

Equation  (137)  is  the  set  of  independent  equations  that  the  non  2ero  eigenvalues  for 
the  control  result  from.  For  the  application  of  pole  placement  algorithms  to  determine  the 
sliding  surfaces,  it  only  has  to  be  recognized  that  it  is  in  the  standard  state  space  format: 

y  =  Ay  +  Bky  O^S) 


with: 

y4  “  ^22 
B  =  -A21 

k  =  sl 


Once  si  is  determined,  the  control  law  can  be  determined  by  substituting  it  into  Equation 
(135). 


Utkin’s  technique  also  allows  for  the  use  of  linear  quadratic  regulator  methods  for 
determination  of  the  sliding  surfaces.  (Utkin  1977)  For  this  method,  it  is  desired  to  minimize  a 
quadratic  performance  index: 


7  =  —  f  x^Qxdt 


(139) 
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where ^  is  a  positive  definite  weighting  matrix.  By  partitioning^  in  the  same  manner  as  A 
was  partitioned  for  equations  (131)and  (132),  Equation  (139)  can  be  rewritten  as: 


/ 


2  J  I-^l  ^2^ Ql\^\  Q\l^l  -^2  222-^2 


(140) 


or: 


{x2^Q*X2 


(141) 


where: 


Q  =  Q22 QixQn  Q\2 
^  -  ^22  ^2\Q\\  Qn 

V  =  jTj  +  !2i/2i2-^i 


With  the  system  on  the  sliding  surfaces,  Equation  (140)  can  be  rewritten  as: 

X2  =  A*  X2+  A21V  (142) 

Equations  (141)  and  (142)  are  in  a  recognizable  form  for  the  application  of  any 
convenient  Unear  quadratic  regulator  solution.  The  Hamiltonian  is: 


(143) 


The  algebraic  Riccati  equation  is: 
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k  +  kA’ -kA^^QuAl^k  +  Q'  =0 


(144) 


The  solution  of  Equation  (144)  results  in: 

V  =  ~Q]  1  A2\kx2 


(145) 


This  result  is  used  with  the  definition  for  v  from  Equation  (141)  to  provide  the  relationship 
between  x\  and  X2.  This  results  in  the  sliding  surface: 

si  =Q^l{Qn+Al,k)  (146) 

With  si  determined,  the  control  law  can  be  determined  by  substitution  into  Equation  (135). 

This  sliding  mode  LQR  controller  design  was  implemented  in  a  MATLAB®  function 
SMLQR.M  which  included  provisions  for  a  QR  factorization  for  the  cases  when  the  B  matrix 
was  not  of  the  form  given  by  Equation  (130).  SMLQR.M  is  included  in  Appendix  A. 

3.  Control  of  chatter 

One  undesirable  aspect  of  sliding  mode  control  is  the  chatter  induced  by  the  nonlinear 
switching  term  near  the  sliding  surfaces.  The  nonlinear  switching  term  in  the  sliding  mode 
control  can  cause  control  chatter  when  the  system  is  near  a  sliding  surface.  One  way  to  reduce 
the  chatter  is  to  use  the  concept  of  a  boundary  layer  around  the  sliding  surface.  First,  a  satsign 
function  is  defined: 


satsignix)  = 


l,x>l 

x-l<x<l 

-l,x<-l 


(147) 


A  boundary  layer  of  thickness  0  around  each  sliding  surface,  is  applied  using  the  safsign 
function: 
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(148) 


/ 

a  I  =  satsign 

\ 


Equation  (147)  can  be  used  to  replace  the  sign  function  used  in  Equation  (129)  with  no  change 
in  the  asymptotic  stability  of  the  system.  There  are  some  effects,  however,  because  the 
dynamics  near  the  sliding  surface  are  not  the  same  as  the  closed  loop  dynamics  which  exist 
when  the  system  is  on  the  sliding  surface.  The  final  control  law  is: 


f 

u  =  Kx+  K^r]satsign  ■ 
\ 


a 

7 


J 


(149) 


Figure  27.  SIMULINK®  model  sliding  mode  controller 


Equation  (149)  was  implemented  as  a  SIMULINK®  model,  shown  in  Figure  27. 


B.  SIMO  SLIDING  MODE  CONTROL  RESPONSE  TO  DISTURBANCES 
When  applied  to  vertical  plane  submarine  control,  sliding  mode  control  has  several 
nuances  which  are  not  obvious  from  inspection  of  the  governing  equations.  In  order  to 
illustrate  these,  the  performance  of  sliding  mode  control  will  be  explored  through  several 
example  cases.  To  keep  the  analytic  derivations  simple,  the  cases  worked  will  be  done  with 
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stern  planes  control  only.  The  general  concepts,  however,  will  be  applicable  to  stern  and  bow 
planes  control. 

The  response  of  sliding  mode  control  to  force  and  moment  disturbances  is 
fundamental  to  its  application  to  submarine  control  in  the  vertical  plane.  These  force  and 
moment  disturbances  could  result  from  a  variety  of  situations.  Examples  of  these  include  out 
of  trim  conditions,  free  surface  effects,  and  wave  forces. 

1,  Basic  sliding  mode  disturbance  response 

The  first  case  will  use  a  basic  sliding  mode  control  law,  that  is  one  without  a 
feedforward  term  or  integral  control.  The  resulting  control  law  is: 

(T  =  5,  w+SjijH- 536  +  54  z  (150) 

o  (151) 

5  =  K]W+  ^2?  +  ^3^  +  'nJ^s^atsign{—) 

(152) 

a )  Linear  analysis  steady  state 

Assuming  that  the  sliding  mode  control  is  not  saturated  and  that  the  control 
deflection  is  less  than  the  maximum,  the  submarine  should  reach  an  equilibrium  state  under  the 
action  of  steady  force  and  or  moment  disturbances.  The  linear  equations  of  motion  in  the 
vertical  plane  with  one  control  are: 


w  =  ax2^q-\- 

(153) 

q  =  a2iWVV  +  a22W^  +  ^23^  +£^2^^^ 

(154) 

e^q 

(155) 

z-w-uO 

(156) 

X  =  wB  +  u 

(157) 

Equations  (150)  through  (157)  can  be  solved  for  the  following  steady  state  condition: 
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(158) 


Fj{-ub^)+  M  j{ub^) 

W^.^.  — - ^ - 2 - 

'+'^2^13  “^1^23 

FA-b,)  +  MAb,) 

— - 2 - 2 

"f“ /72^13  *”^1^23 

_  ^  +ajy  +  u^r}K^(j)b2S,+u^riK,(l)b2S^  +  K^u\ 

U  K^S^(p{b2^\\^  "^^2^13  1^1^21^  b\^2'i^ 

^  K^u'^b^  +auU^  +a,j  +  u^riK^(l)b^Sj  +  K^u^b^  +u^r|K^<j)b^S^ 

^  ^)i‘^40(^2^11^  +^2^13  “^1^21^  ”^1^23) 

g  _  /(j(a2l»^+Q23)+^,2(-Qll“^  -^13) 

'“  +^2^13  “^l‘^2l“^  ~^i%) 


(159) 

(160) 

(161) 


(162) 


If  the  sliding  mode  is  just  saturated,  that  is  |cr  /  ^|  =  1 ,  the  system  will  still  be 
stable,  however  the  gain  parameter  T]  will  be  at  a  critical  value.  If  further  reduced,  the  system 
will  be  unstable.  Assuming  the  sliding  mode  is  just  saturated  this  critical  value,  ,  can  be 
determined. 


K^u^ib^a^^u}  +i>2«i3 -bia^^u^  -^i«23) 
.M  anU^  +a^^  +  Kyb^  +  K,u\ 

K^ib^a^^u^  +b2ai2 -biO^^u^  -^i«23) 


(163) 


For  cases  with  the  gain  parameter  T]  less  than  critical,  a  steady  state  i  results, 
along  with  steady  state  values  in  aU  states  other  than  z  .  By  assuming  a  steady  state  z  and 
solving  Equations  (150)  through  (157)  for  equilibrium,  the  following  equations  result: 


F,, (-a„  - b^uK^ )  +  MdibiuKi  +«,,)  +  u^RMbi^i  1  - &i^i2i ) 

(164) 

K^{b^a2i  -  ^2^11  )  +  “^l(®13^2  “'^23^1) +  ^13^21  '■^23‘*11 

F,I  +b^u^K,-n  +  9„(ai2+a\\U^  +bfU^K^  +b^u^K2) 

(165) 

~  +  aji ) 
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=  u9n  +  z,. 


(166) 


=  0 

-(aii«w,v +fl|30„  +  F^) 

^.v.v  “  ,2 

biU 


(167) 

(168) 


Inspection  of  Equation  (164)  may  cause  the  reader  to  incorrectly  assume  that  a 
nonzero  steady  state  value  of  9  will  exist  in  the  absence  of  disturbances.  This  is  not  the  case, 
however,  because  the  use  of  Equation  (164)  implies  the  lack  of  a  steady  state  z ,  and  therefore 
a  disturbance. 

b)  Nonlinear  analysis  steady  state 

An  analysis  similar  to  that  conducted  on  the  linear  equations  can  be  conducted 
to  determine  the  system  steady  state  response  under  a  constant  disturbance.  The  nonlinear 
equations  of  motion  in  the  vertical  plane  with  one  control  are: 


vi'  =  a,,HW  +  a]2«9  +  ai3  sin(0)  +  l’i«^5  +  cos(0)  +  eii^^  +  ^12^'^ 
q  =  a2iMW  +  U22“9  +^23  sin(^)  +^’2“^^  +  cos(0)  +  e2\q^  +  ^22?^  (^70) 

e  =  q  (171) 

z  =  wcos(0) -Hsin(0)  (172) 

i  =  M'sin(0)  +  t/cos(0)  (173) 


Once  again,  the  basic  sliding  mode  control  of  Equations  (150)  and  (151)  is 
used.  Assuming  that  steady  state  is  possible  (sliding  mode  control  is  not  saturated,  and  the 
control  deflection  is  less  than  the  maximum)  Equations  (169)  through  (173)  subject  to  the  basic 
sliding  mode  control  law  can  be  reduced  to  the  following  nonlinear  equation: 

fli,  sin(0)M^l?2  +^13  sin(0)cos(0)Z72  +  Fj  cos^(6)b2  -  (174) 

Z7,fl2i  sin(0)M^  -^i«23  sin(0)  cos(0)  -  Mj  cos^(6)  -  0 
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Solution  of  Equation  (174)  yields  the  steady  state  value  of  0 .  It  is  of  note  that 
Equation  (1 74)  is  independent  of  the  control  law.  The  correct  root  is  readily  determined  by 
using  the  value  closest  to  the  linear  analysis  (Equation  (160)).  Substitution  of  this  value  into 
Equations  (169)  through  (173)  yields  foUowing  steady  state  results: 


=  «tan(0„) 
4.„=0 


5.  = 


+^13  sin(0„)  +  Fj ) 


by 


(175) 

(176) 

(177) 

(178) 


The  steady  state  value  of  z  given  by  Equation  (178)  is  dependent  upon  the 
control  law  gains.  Following  the  example  of  the  linear  analysis,  a  critical  value  of  the  gain 
parameter,  can  be  determined. 


Ss,  - 

-KAs 

(179) 


Unlike  the  linear  analysis,  the  critical  value  of  the  gain  parameter  is  not  just  a 
linear  combination  of  the  disturbance  forces,  but  rather  requires  nonlinear  solution  for  each 
possible  case. 

If  jet  /  ^i|  >  1  the  sliding  mode  control  will  be  saturated,  and  a  nonzero  steady 
state  z„  will  exist.  For  this  case,  the  control  law  control  law  reduces  to: 


5  =  K^w+  K2q  +  K^O  +  riKygn(a) 


(180) 


Given  the  steady  state  final  condition  in  aU  state  variables  with  the  exception  of 
2,  Equations  (169)  through  (173)  and  (180)  can  be  reduced  to  finding  a  root  of  the  following 
nonlinear  equation: 
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(181) 


{b^K^a2]^ +  biKija2i  — ~^23^1-^l)  I 

b2uK^  +<321 

Fj  cos{9){b2uK^  +ci2\)  ~  cos(0)(fc[MA']  +  fl)])  +  sin(0)(<3[3a2]  —  ^23^11)  _  q 

b2uK\  +<321 


This  can  be  accomplisheti  using  the  using  the  linear  value  of  from  Equation  (1 64)  as  the 


initial  guess.  It  follows  that: 

_  Fj  cos^  (g„ )  +  g,  3  sin(0„ )  cos(g,^ ) 
u{b^uK^+a^^) 

(g;  1  sin(0„ )  + 1>|  kATi  sin(0,, )  + 1>|  cos(0„  )  +  bx  K,t]  cosjO^^ )) 

u(bjuK]  +aii) 


(182) 


w„.  =Mtan(0„)  +  z„ 

(5„  =  A",  w„.  +  ATjO,,  +T]K,sign{o) 


(183) 

(184) 


c)  Disturbance  response  simulation  with  basic  sliding  mode 
These  results  can  then  be  applied  to  the  SUBOFF  hydrodynamic  coefficients. 
For  the  modified  coefficients  used  in  this  thesis,  the  linear  state  space  system  at  six  knots  with 
stem  planes  only  is: 


0.0179 

3.7101 

0.0196 

o' 

W 

■- 0.1 009' 

0.0006 

-0.0680 

-0.0034 

0 

e 

+ 

-.0027 

0 

5  + 

0 

0 

1 

0 

0 

1 

0 

- 10.1269 

0 

z 

0 

0 

(185) 


A  sliding  mode  control  law  is  determined  using  Utkin’s  method.  After  some 
experimentation,  the  diagonal  of  the  minimization  matrix  Q  was  selected  as  Qn  =  100,  Q22  = 
100,  Q33  =  100,  Q44  =  1.  This  yielded  the  following  control  law: 

o  =  \Qw-615591q  - 15.40820  +  0.0830z  (186) 
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(187) 


6  = 


f 

-0.2803W  +  84.8572^  +  7.061 20  +  K.rjsatsign 

V 


\S\  <  0.4 


(188) 


A  moment  disturbance  was  chosen  to  be  controllable  but  give  a  nontrivial 
response.  The  application  of  a  force  of  five  thousand  pounds  and  moment  of  4,573  thousand 
footpounds  resulted  in  a  pure  angular  acceleration  ,  A/^  =0.001  radians/second^.  For  this  state 
space  system,  control  law,  and  disturbance  set,  the  value  was  (0.0661).  Application  of 

the  nonlinear  solutions  yielded  the  steady  state  solutions  in  Table  13. 


T]  I'Hci, 

w^^{Feet  1  Sec) 

q^^iRad  /  sec) 

6, .(Degrees) 

z,,{feet) 

5„  {Radians) 

z^^iFeet  /  Sec) 

0 

-0.9606 

0 

-0.8271 

InSnity 

0.1673 

-0.8143 

0.5 

-1.1791 

0 

-4.3402 

Infinity 

0.1941 

-0.4094 

1.1 

-1.3965 

0 

-7.8518 

-19.5783 

0.2208 

0 

2 

-1.3965 

0 

-7.8518 

-14.6467 

0.2208 

0 

4 

-1.3965 

0 

-7.8518 

-11.6330 

0.2208 

0 

Table  13.  Steady  state  nonlinear  solutions  for  Mj  =0.001  radians/second^ 


The  system  transient  response  was  simulated  using  the  RK45  function  of  the 
SIMULINK®  toolbox.  Figure  28  shows  the  resulting  paths.  The  response  is  given  for  six 
values  of  the  ratio  of  r]  and  .  As  expected,  for  values  of  rj  less  than  critical  the  control 
law  is  unable  to  maintain  a  steady  depth. 
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Figure  28.  Nonlinear  simulation  of  vertical  plane  response  to  a  pure  moment  disturbance 


The  calculations  and  simulations  were  repeated  for  a  pure  force  disturbance.  A 
force  of  43  thousand  pounds  and  a  moment  of  220.4  foot  thousand  pounds  resulted  in  a  pure 
vertical  acceleration  of  0.005  feet/second^.  For  this  state  space  system,  control  law,  and 
disturbance  set,  the  value  ofrj^,,,  was  (0.0466). 


ri  'Tin. 

( Feet  /  Sec 

(Rad  1  sec) 

z„(Feet) 

S  ^^(  Radians) 

z„  (Feet  /  Sec) 

0 

1.5724 

0 

5.4888 

Infinity 

0.2357 

0.5966 

0.5 

1.4156 

0 

2.9632 

Infinity 

0.2549 

0.8902 

1.1 

1.8836 

0 

10.5368 

22.4107 

0.1976 

0 

2 

1-8836 

0 

10.5368 

17.4792 

0.1976 

0 

4 

1.8836 

0 

10.5368 

14.4654 

0.1976 

0 

Table  14.  Steady  state  nonlinear  solutions  for  Fj  =  0.05 
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Figure  29.  Nonlinear  simulation  of  vertical  plane  response  to  a  pure  force  disturbance 

d)  Disturbance  response  with  sliding  mode  feedforward  control 
To  eliminate  the  steady  state  error  induced  by  a  disturbance,  several  techniques 
may  be  employed.  Using  a  feedforward  function  in  the  sliding  mode  control  law  is  one 
technique.  A  feedforward  term  works  by  using  knowledge  of  the  external  disturbance  and 
using  applying  some  degree  of  control  effort.  This  provides  a  steady  state  control  effort  to 
oppose  the  disturbance  without  a  steady  state  error.  This  approach  is  limited  as  it  requires  one 
control  per  zero  error  state  and  may  be  limited  in  other  ways.  A  feedforward  term  can  be 
added  to  the  sliding  mode  control  law  by  changing  the  sliding  surface  to  the  following: 


a  =  S,w+S^q  +  S^e^-S^z  +  S^  (189) 

The  value  of  is  such  that  it  will  equal  the  control  effort  that  is  applied  by  the 
steady  state  quantity  that  is  desired  to  be  zeroed.  Since  the  main  priority  is  to  obtain  zero 
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depth  error,  is  selected  to  equal  the  control  effort  introduced  by  the  previously  calculated 
value  of  steady  state  depth  error: 


5,  =  -tr(5„  -  a:,  w,,  -  K,ej- 
‘  TlK,  ■ 


(190) 


The  steady  state  values  needed  for  Equation  (190)  can  be  determined  from 
linear  or  nonlinear  analysis,  although  the  linear  analysis  will  result  in  a  non-zero  depth  error. 
The  linear  analysis  gives  the  following: 


—  CjMj  +  C2FJ 


(191) 


where: 


C,  =- 


U  J 


^4  J 


Cl^^U  +  ^13^2  “ “”^1^23 


(192) 


C,  = 


a2,+^  +  K^b2  +  K^ub2  +-r^  + 


^2  •  ■  -■■-•■'-j  .^4  j 


fljjW  ^2  ^13^2  “  ^1^21^  ‘”^1^23 


(193) 
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Depth,  Feet 


Figure  31.  Nonlinear  simulation  of  vertical  plane  response  to  a  pure  force  disturbance  with 
a  feedforward  term  based  on  nonlinear  steady  state 

e )  Disturbance  response  with  sliding  mode  integral  control 

Another  means  of  eliminating  steady  state  error  is  by  the  use  of  an  integral 
control  term.  To  accomplish  this,  an  additional  equation  is  added  to  the  state  space 
representation. 


This  forces  a  zero  steady  state  error  in  z  ,  although  there  are  some  additional 
considerations  with  the  use  of  integral  control.  The  resulting  control  law,  with  the  additional 
state,  is: 

CT  =  5, w  +  ^2^  +  ^30  +  54  z  +  S^z,  (1 95) 

CT  (196) 

5  =  K{w  +  +  Ky6  +  K^z  +  riK,satstgn{—) 
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(197) 


|5|  <  0.4 


Based  on  inspection  of  Equations  (1 95)  and  (1 96)  several  conclusions  can  be 
drawn.  First,  because  z  is  included  in  the  proportional  portion  of  the  control  law,  values  of  the 
gain  parameter  which  are  less  than  critical  will  result  in  a  steady  state  error  in  z  for  any 
controllable  disturbance.  Second,  if  there  is  a  steady  state  error  in  z  ,  the  magnitude  of  the 
integral  term,  z, ,  will  tend  to  infinity.  This  can  cause  problems  with  changing  conditions  or 

pathkeeping  as  it  will  delay  the  control  response  to  other  errors. 

At  a  condition  of  steady  state,  with  the  gain  parameter  greater  than  critical,  the 
steady  state  error  in  z  is  zero.  Because  of  this,  the  previously  calculated  values  for  steady  state 
pitch  angle,  heave,  and  control  deflection  are  still  valid  (Equations  (174)  through  (177)). 
Moreover,  because  any  controllable  disturbance  will  result  in  a  steady  state  condition  these 
values  apply  for  cases  where  the  gain  parameter  is  less  than  critical. 

For  conditions  where  the  gain  parameter  is  less  than  critical,  the  resulting 

control  law  is: 

8  =  K^w  + K2q+  K^0+ K^z  +  TjK^  sign{a)  (1 9  8) 

And  the  steady  state  value  of  z  is: 

5,,  -  /iliw,,  -  K2,6,„  -r]K,sign(o)  (199) 


Expressed  in  linear  state  space  form  for  control  law  design,  the  representation 
for  the  SUBOFF  at  six  knots  is: 


w 

■-  0.0179 

3.7101 

0.0196 

0 

O' 

w 

0.1009' 

0.0006 

-0.068 

-  0.0034 

0 

0 

-  0.0027 

M, 

e 

= 

0 

1 

0 

0 

0 

G 

+ 

0 

5  + 

0 

z 

1 

0 

-10.1269 

0 

0 

z 

0 

0 

0 

0 

0 

1 

0 

0 

_  0  _ 

(200) 
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Utkin’s  method  was  applied  to  obtain  a  suitable  control  law.  After  some 
experimentation,  the  diagonal  of  the  minimization  matrix  Q  was  selected  as  Qii  =  100,  Q22  = 
100,  Q33  =  100,  Q44  =  1,  Q55  =  0.01.  This  yielded  the  following  control  law: 


a  =  l.Ow  -  49.8158^  - 12.77250  +  O.lOSOz  +  0.0035z/ 

5  =  -1.6009W+  161.05899  +  24.81250-0.099  lz+ 

|5l  <  0.4 


(201) 

(202) 

(203) 


Simulations  of  the  SUBOFF  under  sliding  mode  integral  control  with  a 
moment  disturbance  are  given  in  Figure  32.  For  this  state  space  system,  control  law,  and 
disturbance  set  the  value  of  t]„,,  was  (0.0484).  Shown  are  five  different  values  the  ratio  of  the 
gain  parameter  to  the  critical  gain  parameter.  For  values  of  the  ratio  larger  than  about  two,  the 
system  exhibited  excessive  oscillation  before  settling  to  zero  depth  error. 

Simulations  of  the  SUBOFF  under  sliding  mode  integral  control  with  a  force 
disturbance  are  given  in  Figure  33.  For  this  state  space  system,  control  law,  and  disturbance  set 
the  value  of  77^,,,  was  (0.0466).  Shown  are  six  different  values  the  ratio  of  the  gain  parameter  to 
the  critical  gain  parameter.  For  values  of  the  ratio  larger  than  about  two,  the  system  exhibited 
excessive  oscillation  before  settling  to  zero  depth  error. 

{)  Sliding  mode  disturbance  response  conclusions 

Submarine  vertical  plane  depth  control  using  sliding  mode  control  can  be 
effectively  achieved  in  the  presence  of  disturbances.  Sliding  mode  control  is  similar  to  linear 
state  feedback  in  that  the  an  external  disturbance  will  result  in  a  steady  state  depth  error. 
However,  if  the  gain  parameters  of  the  sliding  mode  control  are  not  properly  selected,  a  loss  of 
depth  control  can  occur. 

Steady  state  error  can  be  dealt  with  using  feedforward  or  integral  control. 
Integral  control  has  several  advantages.  Application  of  integral  control  does  not  reqioire 
knowledge  of  the  disturbance.  However,  if  the  gain  parameter  is  too  small,  windup  of  the 
integral  term  occurs.  If  the  gain  parameter  is  too  large,  excessive  oscillations  can  be 
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introduced.  The  greatest  advantage  of  integral  control  demonstrated  was  that  gain  parameters 
less  than  critical  did  not  result  in  loss  of  depth  stability. 


X  axis,  Feet 

Figure  32.  Nonlinear  simulation  of  moment  disturbance  using  sliding  mode  integral 

control 

Feedforward  control  exhibited  good  disturbance  compensation,  assuming  that 
the  disturbances  were  measurable.  Given  the  disturbances,  feedforward  values  can  be 
determined  based  upon  the  nonlinear  equations  of  motion,  requiring  periodic  nonlinear  root 
finding,  or  upon  the  linear  solution.  Because  of  the  computational  expense  of  obtaining  the 
nonlinear  solution  and  the  expected  error  in  the  hydrodynamic  coefficients,  a  linear  steady  state 
solution  is  appropriate  for  feedforward  computation. 
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Figure  33.  Nonlinear  simulation  of  force  disturbance  using  sliding  mode  integral  control 


C.  MIMO  SLIDING  MODE  CONTROL  AT  PERISCOPE  DEPTH 

1.  Introduction 

The  purpose  of  using  sliding  mode  control  was  to  provide  an  alternate  means  of 
control  which  relied  upon  all  the  system  states.  The  robust  characteristics  of  sliding  mode 
control  were  thought  to  be  a  good  approximation  to  the  human  operators. 

At  periscope  depth,  experience  dictates  several  desired  conditions.  First,  the  ship  is 
trimmed  heavy  to  counter  the  steady  wave  forces.  Even  more  weight  is  brought  on  after  this 
point  to  allow  for  a  constant  small  positive  trim  angle,  of  several  degrees.  This  provides 
reserve  ballast  which  is  made  available  by  reducing  the  trim  angle.  Finally,  in  sea  state  three,  it 
should  be  very  possible  to  maintain  depth  within  one  foot  of  ordered  depth. 
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The  equations  of  motion  used  for  this  section  are  the  nonlinear  equations  of  submarine 
motion  in  the  vertical  plane.  They  are  different  from  the  equations  used  previously  in  this 
chapter,  as  they  include  both  bow  and  stern  planes.  Also  the  force  and  moment  disturbances 
used  represent  not  only  constant  disturbances,  like  ships  trim,  but  time  varying  wave  forces  as 
well.  Repeated  for  convenience,  the  equations  are: 


w  =  a, ,uw  +  ai2M^  +  ^13  sin(0)  +  Fj  cos(6)  +  (204) 

q  =  a2\uw  +  aiiuq  +<^23  sin(0)  +b2\U^5i,  +^22“^*^.?  cos(0)  +  ^214  +^22^'^  (^0^) 

e  =  g  (206) 

z  =  wcos(0)  -  Msin(0)  (207) 

X  =  wsinC©)  +  ucos(9)  (208) 


2.  Basic  sliding  mode  controller 

The  sliding  mode  controller  is  of  the  same  form  as  before,  although  with  the 
introduction  of  MIMO  control  some  of  the  scalar  terms  become  vectors  or  matrices.  The 
form  of  the  basic  sliding  mode  controller  is: 


CTj  =511^+5129  +  5130  +  514  z 


O2  =  52iW  + 5220  +  ^230 +  ^24  Z 


5^  =  A:iiW+  A'i29+  + 


J 

U.y 


+  ri2K  satsign 


O2 


y<p2  j 


S,  =  K2iW+  K22q  +  K22e  +  rji  K,^^satsign\^- 


+  ri2K,^^satsign 


\h  j 


(209) 

(210) 
(211) 

(212) 

(213) 


(214) 
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As  an  initial  attempt,  Utkin’s  method  was  used  to  determine  a  sliding  mode  control 
law.  After  some  experimentation,  the  diagonal  of  the  minimization  matrix  Q  was  selected  as 
Qii  =  100,  Q22=  100,  Q33  =  100,  Q44  =  1.  This  resulted  in  the  following  control  matrices: 


39.6585  -632.81  -40454  0 

-24.6916  429.65  250.27  0 


(215) 


1 

0 

0 

1 

-00985 

1.7281 

0.0170 

-0.0985 

'10.9711 

-405.06' 

3.083 

252.10 

(216) 


(217) 


The  values  of  the  gain  and  boundary  layer  thickness  parameters  were  taken  as: 


n,  =  m  =  — 
"  ^  20 


(218) 


_  ^2  —  1 


(219) 


A  SIMULINK®  model  was  developed  to  incorporate  the  submarine  dynamics  of 
Chapter  II,  the  wave  forces  of  Chapter  III,  and  the  MIMO  sliding  mode  control  law.  Also 
included  was  the  trim  model  from  Chapter  IV. 

The  model  was  used  to  simulate  a  step  change  in  commanded  depth  from  140  feet  to 
50  feet  in  depth.  To  provide  some  realism  in  the  trim  condition,  the  submarine  was  trimmed 
to  25  thousand  pounds  heavy,  with  no  moment  correction.  Wave  force  values  for  sea  state 
three  were  used,  with  a  relative  heading  of  180  degrees  (head  seas).  Figure  35  shows  the 
resulting  path  taken  by  the  submarine. 
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Figure  34.  SIMULINK®  model  of  submarine  with  wave  forces  and  trim 


At  first  glance,  this  control  scheme  fulfills  most  of  the  desired  characteristics  of  a 
submarine  depth  controller.  Periscope  depth  was  achieved  with  no  overshoot,  and  reasonable 
depth  control  was  maintained.  The  trim  condition  was  selected  so  that  a  steady  state  positive 
trim  angle  would  exist  at  periscope  depth.  During  the  depth  change,  the  maximum  trim  angle 
achieved  was  about  ten  degrees,  which  is  also  very  consistent  with  actual  submarine  practice. 

Inspection  of  Figure  36  shows  some  problems  with  this  particular  controller.  The 
application  of  control  effort  was  excessive.  The  main  reason  for  this  was  the  high  frequency 
variations  in  w  and  ij  induced  by  the  wave  forces.  Because  of  the  combination  of  wave  forces 
and  trim,  the  commanded  depth  was  not  achieved.  The  average  depth  at  periscope  depth  was 
50.75  feet,  as  opposed  to  the  commanded  depth  of  50  feet. 
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Figure  35.  Basic  sliding  mode  performance,  step  change  approach  to  PD 


Although  the  performance  for  the  sea  state  three,  head  seas  was  adequate  with  this 
controller,  it  did  not  perform  well  with  the  other  sea  states  or  headings.  Because  of  this,  it  was 
decided  to  use  this  control  as  a  starting  point  for  a  performance  optimization  for  each  of  the 
four  sea  state  and  direction  cases  available. 

As  was  done  for  state  feedback  control,  the  MATLAB®  CONSTR  function  was  used 
to  perform  the  optimizations.  To  provide  a  general  set  of  design  variables,  the  sliding  mode 
linear  quadratic  regulator  program  was  not  used  to  determine  the  control  law  at  each  step  of 
the  optimization.  Instead,  the  sliding  surface  itself  was  varied  to  change  the  control  law. 
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Figure  36.  State  parameters  for  basic  sliding  mode  approach  to  periscope  depth 


The  formal  optimization  statement  is: 
Minimize: 


> *^32  » *^41  ’*^42  F) 

where: 

z  =  depth  ,  determined  by  nonlinear  simulation 
V 

jiz)dt 

_ 

^mean  ~  ^  J 

V 

H  =  Ballast  added  to  center  of  buojang,  thousands  of  pounds 


(220) 
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F  =  ballast  shifted  from  forward  to  eft,  thousands  of  pounds 


Subject  to: 


real{eigenvalues(A22  -  ^21  ■^2  ))  -  '^max 


(221) 


Deviation  from  the  mean  value  of  depth  vice  the  commanded  was  used  because  of  the 
expected  average  depth  error. 

This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 


I _ xfeet _ xfeet _ | 

Figure  37.  Simulation  with  basic  sliding  mode  control  in  sea  state  three  (head  sea  direction) 
seas),  the  optimized  response  is  shown  in  Figure  37.  The  results  of  the  four  optimizations  are 

shown  in  Table  15.  For  the  RMS  error  and  maximum  error,  the  optimized  values  are  given, 

along  with  their  percentage  of  the  initial  values.  In  aU  cases,  use  of  the  optimization  resulted  in 

reduction  of  the  mean  square  depth  error  (measured  from  the  average  depth).  Reduction  of 

the  maximum  error  was  also  achieved. 
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Sea  State /Direction 

3 /head 

3/beam 

4/head 

4/ beam 

Initial  Values 

s 

1 

0 

1 

0 

1 

0 

1 

0 

0 

1 

1 

0 

1 

0 

1 

-0.0985 

1.7281 

1.7281 

-0.1388 

1.377 

-0.0932 

1.7281 

0.0170 

-0.0985 

0.0170 

-0.0985 

0.0191 

-0.694 

0.0170 

-0.0903 

KT 

39.658 

-24.69 

39.658 

-24.69 

27.88 

-17.34 

36.32 

-22.61 

-632.81 

429.6 

-632.81 

429.6 

-491.0 

341.0 

-632.77 

429.7 

-404.54 

250.27 

-404.54 

250.27 

-285.23 

175.8 

-370.71 

229.21 

0 

0 

0 

0 

0 

0 

0 

0 

m  /J?2 

0.05/0.05 

0.05/0.05 

0.05/0.05 

0.0445/0.0445 

H/F  (lO-"^  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

53.46 

49.55 

52.83 

54.53 

RMS  Error  (feet) 

1.48 

3.37 

1.63 

2.75 

Maximum  Error  (feet) 

3.00 

8.03 

4.20 

9.78 

Sliding  surface 

-0.8725  +  0.50631 

-0.8725  +  0.50631 

-0.6981  +  0.48201 

-0.8726  +  0.41721 

eigenvalues 

-0.8725 

-  0.5063i 

-0.8725 

-  0.5063i 

-0.6981 

-  0.4820i 

-0.8726 

-  0.41 72i 

Optimized  Values 

SMii 

S 

1 

0 

1 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

-0.0816 

1.7244 

-0.0196 

1.884 

-0.1525 

1.3761 

-0.933 

1.7281 

0.0205 

-0.1203 

0.0037 

-0.0208 

0.0176 

-0.0759 

0.0170 

-0.0903 

¥J 

48.50 

-30.16 

8.008 

-5.124 

30.475 

-18.97 

36.322 

-22.615 

-631.1 

428.78 

-695.15 

469.2 

-490.8 

340.74 

-632.76 

429.67 

-494.0 

305.61 

-84.01 

52.11 

-311.5 

192.3 

-370.75 

229.24 

0 

0 

0 

0 

0 

0 

0 

0 

Hi  /I?2 

0.0488/0.0510 

0.0457/0.0457 

0.0501/0.0501 

0.0446/0.0445 

H/F  (103  pounds) 

19.9/0 

19.6/0.0 

20.0/0 

20.0/0 

Mean  Depth  (feet) 

55.44 

57.61 

55.53 

54.57 

RMS  Error  (feet) 

0.27  (18%) 

0.84  (25%) 

0.683  (42%) 

1.43  (52%) 

Maximum  Error  (feet) 

0.77  (26%) 

2.36  (29%) 

1.92  (46%) 

4.14  (42%) 

Sliding  surface 

-0.8725  +  0.69451 

-1.7649 

-0.6968  +  0.54341 

-0.8726  +  0.41731 

eigenvalues 

-0.8725 

-  0.6945i 

-0.1229 

-0.6968 

-  0.5434i 

-0.8726 

-  0.41731 

Table  15.  Optimized  basic  sliding  mode  control  law  results  and  performance 
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3.  Disturbance  feedforward 

The  sliding  mode  control  can  be  implemented  with  a  disturbance  feedforward  to 
correct  average  depth  error.  This  can  be  implemented  inside  or  outside  of  the  sliding  surface. 
For  this  example,  the  disturbance  feedforward  was  implemented  inside  the  sliding  surface 
calculation. 


(7  =  5  jc  +  iSc 


(222) 


Assuming  that  neither  sliding  surface  is  saturated,  that  control  deflection  is  within 
limits,  and  using  a  linear  analysis,  the  steady  state  value  of  the  depth  error  can  be  written  as  a 
linear  combination  of  the  force  and  moment  disturbances  (Appendix  B). 


^ss  ^commanded  ^i^d'^^2^d 


(223) 


To  eliminate  the  depth  error,  it  is  necessary  to  apply  the  same  amount  of  control  effort  that  the 
steady  state  error  provides  within  the  sliding  surface: 

■54.1  (224) 


^5  ^commanded ) 


■’42 


This  results  in  the  following  control  law: 


K,2 

^13 

0‘ 

^  ^  commanded 

/21 

K22 

^23 

0 

^ 

^  ^commanded 

^  ~  ''^commanded 


^  ^  commanded 


K 


•'ll 


CJ2 


'  1 

0  ■ 

T 

^  ^commanded 

0 

1 

^  commanded 

■^31 

*^32 

^ 

^  ^commanded 

.•=^41 

‘^42. 

_  ^ ^commanded 

K.. 


where  the  force  and  moment  disturbances  are  filtered. 


rj^satsigm 


‘t>\ 


.  f  0-2 

T]2satsign 

\V2  j 


^\^4\ 

^2^4\ 

^2^42^ 

(225) 


(226) 


(227) 
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The  initial  sliding  surface  and  gains  from  the  basic  sliding  mode  control  law  was  used 
as  the  starting  point  for  optimization.  The  formal  optimization  statement  is: 


^(*^31  >  ■^32  » *^41  » *^42  ’  ^co  >  >  ^2  »  “ 


j 

1  ”  ^commanded  ) 


(228) 


Subject  to: 


real{eigenvalues{A22  -  ^21 ))  - 


(229) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimized  response  is  shown  in  Figure  38.  The  results  of  the  four  optimizations  are 
shown  in  Table  16.  For  the  RMS  error  and  maximum  error,  the  optimized  values  are  given, 
along  with  their  percentage  of  the  initial  values. 
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Sea  State /Direction 
Initial  Values 


3/head 


3 /beam 


4/head 


4/beam 


Initial  Values 

S'!?!! 

S 

1 

0 

1 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

1.7281 

-0.0985 

1.7281 

1.7281 

-0.0985 

1.7281 

0.0170 

-0.0985 

0.0170 

0.0170 

-0.0985 

0.0170 

-0.0985 

KJ 

39.658 

-24.69 

39.658 

-24.69 

39.658 

-24.69 

39.658 

-24.69 

-632.81 

429.6 

-632.81 

429.6 

-632.81 

429.6 

-632.81 

429.6 

-404.54 

250.27 

-404.54 

250.27 

-404.54 

250.27 

-404.54 

250.27 

0 

0 

0 

0 

0 

0 

0 

0 

0)^^,  (radians/ second) 

0.25 

0.25 

0.25 

0.25 

m  fvi 

0.05/0.05 

0.05/0.05 

0.05/0.05 

0.05/0.05 

H/F  (10^  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

55.28 

54.99 

55.10 

53.41 

RMS  Error  (feet) 

0.4181 

1.50 

1.01 

7.04 

Maximum  Error  (feet) 

1.225 

3.85 

3.50 

18.75 

Sliding  surface 

-0.8725  +  0.5063i 

-0.8725  +  0.5063i 

-0.8725  +  0.5063i 

-0.8725  +  0.5063i 

eigenvalues 

-0.8725 

-  0.5063i 

-0.8725 

-  0.5063i 

-0.8725 

-  0.5063i 

-0.8725 

-  0.5063i 

Optimized  Values 

S 

1 

0 

1 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

-0.0981 

1.7281 

-0.1255 

2.1291 

3.8086 

-0.1167 

1.8775 

0.0166 

0.096 

0.0400 

-0.0473 

-0.1476 

-0.0240 

KJ 

38.67 

-24.08 

19.16 

-11.70 

58.9 

-37.3 

9.906 

-5.992 

-632.8 

429.65 

-795.53 

530.65 

-1496.0 

948.4 

-693.54 

467.3 

-394.5 

244.1 

-196.91 

118.74 

-599.1 

377.5 

-95.02 

60.9 

0 

0 

0 

0 

0 

0 

0 

0 

(radians/ second) 

0.25 

0.350 

0.267 

0.245 

m  frii 

0.05/0.05 

0.0542/0.0578 

0.1853/0.0544 

0.0596/0.0548 

H/F  (10^  pounds) 

20/0 

20.2/-0.1 

20.7/1.3 

20.1/0.0 

Mean  Depth  (feet) 

55.22 

55.47 

55.45 

2.00 

RMS  Error  (feet) 

0.405  (97%) 

1.22  (81%) 

1.01  (100%) 

1.99  (28%) 

Maximum  Error  (feet) 

1.36  (111%) 

3.20  (83%) 

2.52  (72%) 

4.93  (26%) 

Sliding  surface 

-0.8723  +  0.4813i 

-1.8706 

-3.4649 

-1.7409 

eigenvalues 

-0.8723 

-  0.481 3i 

-0.2984 

-0.2992 

-0.1211 

Table  16.  Optimized  sliding  mode  control  with  disturbance  feedforward  results  and 


performance 
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4.  Integral  control 

The  sliding  mode  control  law  can  be  augmented  with  integral  control  on  depth  to 
remove  the  average  depth  error.  This  results  in  the  following  control  law: 


^bp 

'Kn 

Kx2 

0] 

KJ 

_^21 

^22 

^23 

K2, 

oj 

^  ^commanded 
^  Q  commanded 
^  ^commanded 
^  ~  ^commanded 

Zi 


K.. 


'h] 


K 


hi 


r  ^  \ 


7],  satsigri^ 
r}2satsign\ 


r 


02 


(230) 


’  1 

0  ’ 

T 

^commanded 

0 

1 

^  ^  commanded 

*^31 

•^32 

^  ^commanded 

*^42 

^  ^commanded 

_*^51 

*^52. 

|5|<5, 


(231) 


(232) 


After  a  stable  set  of  gains  was  determined,  the  controller  was  optimi2ed  to  niinimi2e 
the  deviation  from  the  commanded  depth.  The  formal  optimi2ation  statement  is: 


Minimi2e: 


^(■^3 1  > -^32  >  “^41  >  *^42  >  •^51  > '^52  >  ^1  >  ^2  •  ^>  ^) 


I 


J 

1  ^commanded  ) 


(233) 


Subject  to: 


real{eigenvalues{A22  —  ))  - 


(234) 


This  approach  was  used  for  each  of  the  four  sea  state  cases.  For  sea  state  three  (head 
seas),  the  optimi2ed  response  is  shown  in  Figure  39.  The  results  of  the  four  optimi2ations  are 
shown  in  Table  17. 
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Sea  State/Direction 


Initial  Values 


3 /head 


3/beam 


4/head 


s 

1 

0 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

1 

-0.1038 

1.7582 

-0.1038 

1.7582 

-0.1038 

1.7582 

-0.1038 

1.7582 

0.0179 

-0.1038 

0.0179 

-0.1038 

-0.1038 

0.0179 

-0.1038 

0.0005 

-0.0031 

0.0005 

-0.0031 

-0.0031 

0.0005 

-0.0031 

KT 

41.79 

-26.01 

41.79 

-26.01 

41.79 

-26.01 

41.79 

-26.01 

-645.08 

437.2 

-645.08 

437.2 

-645.08 

437.2 

-645.08 

437.2 

-426.14 

263.63 

-426.14 

263.63 

-426.14 

263.63 

-426.14 

263.63 

1.268 

-0.7835 

1.268 

-0.7835 

1.268 

-0.7835 

1.268 

-0.7835 

0 

0 

0 

0 

0 

0 

0 

0 

m  /ri2 

0.05/0.05 

0.05/0.05 

0.025/0.025 

0.05/0.05 

H/F  (10’  pounds) 

20/0 

20/0 

20/0 

20/0 

Mean  Depth  (feet) 

54.82 

54.80 

55.00 

54.73 

RMS  Error  (feet) 

0.4671 

1.47 

0.707 

13.3 

Maximum  Error  (feet) 

1.14 

4.23 

2.03 

31.62 

Sliding  surface 

-0.8722  +  0.5063i 

-0.8722  +  0.5063i 

-0.8722  +  0.5063i 

-0.8722  +  0.5063i 

eigenvalues 

-0.8722 

-  0.5063i 

-0.8722  -  0.5063i 

-0.8722  -  0.5063i 

-0.8722  -  0.5063i 

-0.0316 

-0.0316 

-0.0316 

-0.0316 

Optimized  Values 


S 


7.931 

-32.52 

-83.24 

0.3965 

0 


0 

1 

0.2485 

-0.0202 

-0.0010 


-4.920 

56.91 

50.04 

-0.245 

0 


30.73 

■580.47 

314.13 

1.266 

0 


0 

1 

1.5995 

-0.0766 

-0.0031 


-19.18 

397.3 

194.52 

-0.783 

0 


1 

0 

-0.1428 

0.0177 

0.0003 


25.0 

-1469.6 

-256.2 

0.20 

0 


0 

1 

3.7927 

-0.0624 

0.0006 


-15.6 

950.0 

157.8 

0.10 

0 


/^2 


H/F  (105  pounds) 


Mean  Depth  (feet) 


RMS  Error  (feet) 


Maximum  Error  (feet) 


Sliding  surface 
eigenvalues 


0.0468/0.0464 


19.6/-0.1 


54.98 


0.2345  (50%) 


0-789  (69%) 


-0.1090  +  0.4314i 
-0-1090 -0.4314i 
-0.0500 


0.0461/0.0453 


19.6/0.0 


55.01 


0.713  (49%) 


1.90  (45%) 


-0.2418  +  0.28711 
-0.2418-0.28711 
-0.0481 


0.025/0.025 


20.0/0.0 


55.02 


0.693  (98%) 


1.92  (95%) 


-0.7834  +  0.3260i 
-0.7834  -  0.3260i 
-0.0446 


0.0436/0.0553 


19.4/-0.1 


55.16 


1.47  (11%) 


4.15  (13%) 


-3.6194 

-0.1971 

0.0061 


Table  17.  Optimized  sliding  mode  integral  control  law  results  and  performance 
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D.  CONCLUDING  REMARKS 

A  comparison  between  the  quality  of  control  achieved  by  sliding  mode  control  could 
conclude  that  the  sliding  mode  control  was  inferior  to  state  feedback  control.  This  comparison 
would,  however,  neglect  the  added  benefits  of  sliding  mode  control.  The  robust  character  of 
sliding  mode  control  with  the  ability  to  provide  reliable  control  for  submarine  control  given 
uncertain  hydrodynamic  coefficients  has  been  demonstrated  for  the  NPS  autonomous 
underwater  vehicle  program  (Hawkinson,  1990). 

The  sliding  mode  optimizations  did  not  substantially  reduce  the  control  chatter  and 
attendant  high  actuation  rates.  Variations  of  the  sliding  mode  boundary  thickness  did  not 
alleviate  the  chatter. 

Table  18  gives  a  summary  of  the  RMS  error  achieved  by  each  of  the  sliding  mode 
optimizations.  For  comparison,  it  also  includes  the  full  state  feedback  results  from  Chapter  TV. 
Although  these  were  larger  than  the  corresponding  full  state  feedback  cases,  the  sliding  mode 
control  proved  to  be  much  more  robust  in  response  to  step  changes  in  commanded  depth.  The 
sliding  surface  eigenvalues  exhibited  much  more  damping  than  the  corresponding  cases  of  full 
state  feedback  control.  Also,  it  seemed  to  provide  a  more  realistic  average  pitch  angle  for 
periscope  depth  operations. 


Sea  State/Direction 

3/Head 

3/Beam 

4/Head 

4/ Beam 

Control  Scheme 

Full  State 

0.037 

0.2638 

0.2683 

1.24 

Basic  sliding  mode 

0.27 

0.84 

0.683 

1.43 

Full  State  with  feedforward 

0.0928 

0.4121 

0.400 

0.792 

Sliding  mode  with  disturbance  feedforward 

0.405 

1.22 

1.01 

1.99 

Full  State  integral 

0.0414 

0.372 

0.536 

1.96 

Sliding  mode  with  integral  control 

0.2345 

0.713 

0.693 

1.47 

Table  18.  Optimized  RMS  error  (feet)  of  sliding  mode  control  and  full  state  feedback 
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VI.  GRAPHICAL  DISPLAY 


A.  INTRODUCTION 

In  conducting  ship  control  at  periscope  depth,  a  submarine  diving  officer  relies  on  a 
variety  of  indications,  meters,  and  some  verbal  reports  to  maintain  the  ship  within  the  required 
depth  band.  In  addition  to  the  displayed  parameters,  the  ship’s  control  party  also  has  their 
inertial  reference,  or  “the  seat  of  the  pants”.  It  is  perhaps  the  inertial  reference  which 
differentiates  between  great  ship’s  control  parties,  and  the  merely  adequate. 

A  submarine  diving  officer  must  track  status  of  many  ship’s  systems  in  addition  to  ship 
control.  The  items  which  the  diving  officer  must  monitor  include: 

•Mast  positions 

•Proximity  of  any  portion  of  the  ship  to  broaching  (sail,  rudder,  mast  fairing) 

•Water  depth  (general  terms) 

•Ship’s  relationship  to  the  submerged  operating  envelope 
•Trim 

•Speed 

•Water  density 

•Ship’s  evolutions  (trash  disposal,  ventilation,  etc.) 

•Towed  array,  floating  wire  antenna 

In  many  cases,  the  tracking  tool  most  used  is  the  diving  officer’s  mental  picture. 
Unfortunately,  the  ability  to  keep  a  clear  status  on  many  issues  varies  with  fatigue  and  among 
individuals.  This  chapter  gives  the  current  conditions  of  the  interface  between  the  diving 
officer  and  ship’s  control,  and  proposes  a  different  display  medium  to  improve  operations. 

B.  CURRENT  DIVING  OFFICER  INTERFACE 

To  maintain  a  complete  status,  the  diving  officer  has  few  tools  at  his  disposal.  He  must 
rely  on  looking  around  at  several  different  panels  to  get  mast  status,  soundings,  and  water 
density  while  supervising  the  planesman.  If  an  unplanned  event,  for  example  broaching, 
occurs  the  only  record  for  reconstruction  is  the  memory  of  the  operators. 
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'Ihe  gauges  and  meters  used  for  ship  control,  while  designed  with  generally  appropriate 
time  constants,  are  not  adaptable  for  a  given  circumstance.  For  the  most  part,  the  same 
indications  are  used  for  high  speed  transit,  periscope  depth,  and  tactical  operations.  Figure  40 
and  Figure  41  show  some  of  the  indications  used  on  board  the  USS  Nautilus  (SSN  571). 
Although  Nautilus  is  now  a  museum,  the  design  of  submarine  ship’s  control  panels  has  not 
changed  significantly. 


Figure  40.  USS  Nautilus  planes  position  indications 


Figure  41.  USS  Nautilus  pitch  angle  indication 


The  current  sys  tem  of  ship  status  display  is  very  reliable,  with  redundant  indications  for 
important  items  and  some  purely  mechanical  indicators.  It  falls  short  in  the  area  of  presenting  an 
integrated  status.  It  is  the  writers  opinion  that  the  display  system  degrades  the  performance  of 
the  ship’s  control  party.  With  some  parameters  not  displayed  and  others  not  conditioned  for  the 
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ship’s  operating  mode,  operation  near  the  surface  in  heavy  sea  states  is  extremely  difficult. 

Skilled  operators  rely  upon  the  existing  indications,  as  well  as  “the  seat  of  their  pants”  to 
maintain  depth.  Even  so,  it  is  a  solid  accomplishment  to  keep  from  broaching  during  extended 
periscope  depth  operations. 

Even  more  complex  are  operations  in  shallow  water,  llie  proximity  to  grounding 
complicates  all  aspects  of  ship  control.  The  diving  officer  must  be  constantly  aware  of  the  water 
beneath  the  keel  available  for  casualty  recovery.  Because  nonzero  pitch  angles  will  cause  one  end 
of  the  ship  to  be  deeper  than  indicated,  this  must  also  be  accounted  for. 

C.  PROPOSED  DISPLAY 

To  incorporate  the  desired  indications  in  a  single  display,  a  radically  different  approach  is 
taken.  Rather  than  rely  on  meters  and  gauges  for  the  state  of  the  ship,  a  screen  is  used.  Figure 
42  shows  the  proposed  display.  A  crude  version  of  this  display  was  developed  using  the 
SIMULINK®  Animation  Toolbox®  (Figure  43). 


Towed  Array 


Last  Sounding 


Figure  42.  Proposed  graphical  display  of  submarine  control  status 
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Figure  43.  SIMULINK®  animation  of  depth,  pitch  angle,  and  planes  angles 


By  integrating  the  ship  status  into  one  display,  numerous  improvements  can  be 
realized.  The  relationship  of  the  submarine  to  the  bottom  and  the  surface  is  clearly  shown. 
With  the  bottom  contour  information  from  a  database,  the  diving  officer  has  a  continuous 
sense  of  the  ship’s  proximity  to  grounding.  In  addition  as  sounding  data  is  obtained,  it  can  be 
displayed. 

The  use  of  a  digital  display  paradigm  allows  the  display  to  be  modified  to  support 
different  operating  modes.  Because  of  the  relationship  between  ship  control  and  safety,  the 
settings  would  be  chosen  based  on  a  commanding  officer  approved  doctrine.  This  would 
allow  the  operators  to  adjust  the  display  system  to  best  fit  needs,  and  adapt  it  to  new 
circumstances  or  missions.  Alerts  and  alarms  could  readily  added  as  the  situation  warranted. 

To  assist  the  diving  officer  in  maintaining  status  on  the  wave  forces,  several  bar  graphs 
were  added  to  show  the  net  force  that  the  ship’s  angle  and  planes  were  applying  at  a  given 
time.  These  quantities  would  be  filtered  to  provide  a  relevant  average.  Provided  the  averaging 
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interval  was  appropriate,  this  would  queue  the  diving  officer  to  order  trim  changes  in  response 
to  changing  environmental  conditions. 

The  periscope  video  in  the  bottom  left  hand  corner  would  provide  critical  feedback  for 
the  dive.  This  would  make  the  scope’s  position  relative  to  the  surface  apparent  (another 
indication  of  depth),  and  allow  the  diving  officer  to  be  somewhat  aware  of  the  tactical 
situation.  A  close  or  new  contact  would  prompt  the  diving  officer  to  review  mast  exposure, 
which  is  also  on  the  same  display. 

Safety  of  shallow  water  operations  would  be  enhanced  by  presenting  a  clear  picture  of 
the  ship’s  relationship  to  the  bottom.  During  evasive  action,  the  ship’s  control  party  and  the 
Officer  of  the  Deck  would  be  working  with  common  knowledge  of  available  water  beneath  the 
keel,  and  the  contour  ahead  of  the  ship. 

Ship’s  status  could  be  recorded,  to  allow  playback  for  the  reconstruction  of  unplanned 
events.  Figure  44  shows  a  possible  data  architecture  to  support  the  display. 
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Periscope 

video 


Figure  44.  Graphical  display  data  paths 

D.  CONCLUDING  REMARKS 

The  integration  of  the  pertinent  ship  parameters  in  one  display  should  yield  dramatic 
improvements  in  submarine  periscope  depth  operations.  The  diving  officers  improved 
awareness  should  reduce  fatigue  levels,  allow  for  slightly  lower  speeds  for  a  given  sea  state 
(reducing  mast  feather),  and  enable  a  much  more  complete  environmental  picture  for  the  ship’s 
control  party.  This  awareness  should  increase  the  confidence  of  the  ship’s  control  party  during 
demanding  shallow  water  operations,  reduce  the  likelihood  of  grounding  or  broaching,  and 
provide  an  improved  level  of  support  for  the  Officer  of  the  Deck. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

As  additional  states  were  added  to  the  control  laws,  the  level  achieved  by  the  control 
law  optimizations  generally  improved.  For  the  full  state  feedback  control  cases,  the  depth 
control  exceeded,  in  the  authors  experience,  what  is  achievable  by  manual  control. 

Ship’s  control  parties  use  more  information  about  the  state  of  the  ship  than  is  available 
from  the  explicit  indications.  Also,  the  success  of  the  disturbance  feedforward  control  suggests 
that  averaged  net  force  and  moment  would  be  of  value  to  the  ship’s  control  party. 

Sliding  mode  control  provided  well  damped  dynamics,  with  a  robust  behavior  in 
response  to  changes  in  commanded  depth.  Although  sliding  mode  control  did  not  achieve  the 
very  low  RMS  errors  of  full  state  feedback,  the  gains  which  it  employed  were  smaller  and  more 
realistic. 

B.  RECOMMENDATIONS 

Although  very  good  depthkeeping  was  achieved  with  state  feedback  control,  the 
optimization  schemes  used  resulted  in  control  surface  chatter,  and  very  high  control  surface 
rates.  To  improve  the  quality  of  the  model  and  provide  more  realistic  planesman  action 
several  features  could  be  added  to  the  control  laws  and  optimization  routines: 

•  Investigate  other  sea  states,  speeds  and  sea  directions 

•  Incorporate  control  surface  rate  limits 

•  Include  control  surface  chatter  in  the  optimization  objective  functions 

•  Use  of  Kalman  filtering  to  provide  state  estimation  and  filtering 

•  Investigate  the  use  of  depth  rate  for  feedback  control  in  place  of  heave 

The  application  of  a  new  display  system  to  an  operating  submarine  is  a  major 
undertaking.  Recommended  steps  to  find  a  new  manual  submarine  depth  control  paradigm 
are: 


•  Application  of  system  identification  techniques  to  submarine  operating  data  to 
investigate  the  nature  of  the  human  control 
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•  Trials  of  a  display  onboard  an  appropriate  submarine  and  or  a  submarine  dive  trainer 

•  Use  of  recorded  submarine  operating  data  to  provide  for  “instant  replay”  training  of 
ship’s  control  personnel 
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APPENDIX  A 


All  computer  code  and  SIMULINK®  models  used  in  this  thesis  are  available  from 
Professor  Fotis  PapouJias,  Naval  Postgraduate  School.  The  computer  programs  and 
SIMULINK®  models  used  were: 

Programs: 

SUBOFF.M'-  initializes  SUBOFF  hydrodynamic  coefficients 
AXNL.M’-  performs  nonlinear  Ax+E[qq;wq]  calculation 

SMLQRbM-  determines  MIMO  sliding  mode  control  law  using  Utkin’s  method  (LQR) 

WF_INI.M-  reads  wave  data  files,  processes 

WFORCE.M-  calculates  wave  forces  for  a  given  depth  and  time 

SB_INLM-  initializes  model  variables  for  MIMO  vertical  plane  submarine  model 

SBI_INLM-  initializes  model  variables  for  MIMO  vertical  plane  submarine  model  with  integral 

depth  control 

SB_SM.M-  calculates  MIMO  SM  control  law  from  Q  matrix 
SB_SS.M-  calculates  MIMO  SM  control  law  from  sliding  surface 
SB_SMFF.M-  determines  the  feed  forward  matrix  for  a  given  sliding  mode  control  law 
SB_PD.M-state  feedback  control  law 

SB_PDFF.M-  determines  the  feed  forward  matrix  for  a  state  feedback  control  law 
OBJ2.M  -  Objective  function  for  pitch  /  depth  feedback  optimizations 
OPT2A.M-  Optimization  program  for  OBJ2.M  and  sea  state  three  (head  seas) 

OPT2B.M-  Optimization  program  for  OBJ2.M  and  sea  state  three  (beam  seas) 

OPT2C.M-  Optimization  program  for  OBJ2.M  and  sea  state  four  (head  seas) 

OPT2D.M-  Optimization  program  for  OBJ2.M  and  sea  state  four  (beam  seas) 

OBJ2ff.M  -  Objective  function  for  pitch  /  depth  feedback  with  disturbance  feedforward 
optimizations 

OPT2FFA.M-  Optimization  program  for  OBJ2FF.M  and  sea  state  three  (head  seas) 
OPT2FFB.M-  Optimization  program  for  OBJ2FF.M  and  sea  state  three  (beam  seas) 
OPT2FFC.M-  Optimization  program  for  OBJ2FF.M  and  sea  state  four  (head  seas) 


’  Given  after  list  of  programs 


115 


OPT2FFD.M-  Optimization  program  for  OBJ2FF.M  and  sea  state  four  (beam  seas) 

OBJ2I.M  -  Objective  function  for  pitch  /  depth  feedback  with  integral  depth  control 
optimizations 

OPT2IA.M-  Optimization  program  for  OBJ2LM  and  sea  state  three  (head  seas) 

OPT2IB.M-  Optimization  program  for  OBJ2I.M  and  sea  state  three  (beam  seas) 

OPT2IC.M-  Optimization  program  for  OBJ 21. M  and  sea  state  four  (head  seas) 

OPT2ID.M-  Optimization  program  for  OBJ2I.M  and  sea  state  four  (beam  seas) 

OBJ3.M  -  Objective  function  for  full  state  partial  distribution  feedback  optimizations 
OPT3A.M-  Optimization  program  for  OBJ3.M  and  sea  state  three  (head  seas) 

OPT3B.M-  Optimization  program  for  OBJ3.M  and  sea  state  three  (beam  seas) 

OPT3C.M-  Optimization  program  for  OBJ3.M  and  sea  state  four  (head  seas) 

OPT3D.M-  Optimization  program  for  OBJ3.M  and  sea  state  four  (beam  seas) 

OBJ3FF.M  -  Objective  function  for  full  state  feedback  with  disturbance  feedforward 
optimizations 

OPT3FFA.M-  Optimization  program  for  OBJ3FF.M  and  sea  state  three  (head  seas) 
OPT3FFB.M-  Optimization  program  for  OBJ3FF.M  and  sea  state  three  (beam  seas) 
OPT3FFC.M-  Optimization  program  for  OBJ3FF.M  and  sea  state  four  (head  seas) 
OPT3FFD.M-  Optimization  program  for  OBJ3FF.M  and  sea  state  four  (beam  seas) 

OBJ3LM  -  Objective  function  for  full  state  feedback  with  integral  depth  control  optimizations 
OPT3IA.M-  Optimization  program  for  OBJ3LM  and  sea  state  three  (head  seas) 

OPT3IB.M-  Optimization  program  for  OBJ3I.M  and  sea  state  three  (beam  seas) 

OPT3IC.M-  Optimization  program  for  OBJ3I.M  and  sea  state  four  (head  seas) 

OPT3ID.M-  Optimization  program  for  OBJ3LM  and  sea  state  four  (beam  seas) 

OBJ3.M  -  Objective  function  for  full  state  feedback  optimizations 
OPT4A.M-  Optimization  program  for  OBJ4.M  and  sea  state  three  (head  seas) 

OPT4B.M-  Optimization  program  for  OBJ4.M  and  sea  state  three  (beam  seas) 

OPT4C.M-  Optimization  program  for  OBJ4.M  and  sea  state  four  (head  seas) 

OPT4D.M-  Optimization  program  for  OBJ4.M  and  sea  state  four  (beam  seas) 

OBJ4FF.M  -  Objective  function  for  full  state  partial  distribution  feedback  with  disturbance 
feedforward  optimizations 

OPT4FFA.M-  Optimization  program  for  OBJ4FF.M  and  sea  state  three  (head  seas) 
OPT4FFB.M-  Optimization  program  for  OBJ4FF.M  and  sea  state  three  (beam  seas) 
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OPT4FFC.M-  Optimization  program  for  OBJ4FF.M  and  sea  state  four  (head  seas) 
OPT4FFD.M-  Optimization  program  for  OBJ4FF.M  and  sea  state  four  (beam  seas) 

OBJ4LM  -  Objective  function  for  full  state  partial  distribution  feedback  with  integral  depth 
control  optimizations 

OPT4IA.M-  Optimization  program  for  OBJ4LM  and  sea  state  three  (head  seas) 

OPT4IB.M-  Optimization  program  for  OBJ4LM  and  sea  state  three  (beam  seas) 

OPT4IC.M-  Optimization  program  for  OBJ4I.M  and  sea  state  four  (head  seas) 

OPT4ID.M-  Optimization  program  for  OBJ4LM  and  sea  state  four  (beam  seas) 

OBJ7.M  -  Objective  function  for  sliding  mode  control  optimizations 
OPT7A.M-  Optimization  program  for  OBJ7.M  and  sea  state  three  (head  seas) 

OPT7B.M-  Optimization  program  for  OBJ7.M  and  sea  state  three  (beam  seas) 

OPT7C.M-  Optimization  program  for  OBJ7.M  and  sea  state  four  (head  seas) 

OPT7D.M-  Optimization  program  for  OBJ7.M  and  sea  state  four  (beam  seas) 

OBJ7FF.M  -  Objective  function  for  sliding  mode  control  with  disturbance  feedforward 
optimizations 

OPT7FFA.M-  Optimization  program  for  OBJ7FF.M  and  sea  state  three  (head  seas) 
OPT7FFB.M-  Optimization  program  for  OBJ7FF.M  and  sea  state  three  (beam  seas) 
OPT7FFC.M-  Optimization  program  for  OBJ7FF.M  and  sea  state  four  (head  seas) 
OPT7FFD.M-  Optimization  program  for  OBJ7FF.M  and  sea  state  four  (beam  seas) 

OBJ7I.M  -  Objective  function  for  sliding  mode  control  with  integral  depth  control 
optimizations 

OPT7IA.M-  Optimization  program  for  OBJ7LM  and  sea  state  three  (head  seas) 

OPT7IB.M-  Optimization  program  for  OBJ7LM  and  sea  state  three  (beam  seas) 

OFT/IC-M-  Optimization  program  for  OBJ7I.M  and  sea  state  four  (bead  seas) 

OPT7ID.M-  Optimization  program  for  OBJ7I.M  and  sea  state  four  (beam  seas) 

OBJ7.M  -  Objective  function  for  full  state  feedback  optimizations 

Models; 

SMSW-M-  MIMO  sliding  mode  control  submarine  control  model,  with  wave  forces  and  trim 

SMSWFF.M-  Same  as  SMSW.M,  with  disturbance  feedforward 

SMT.M-  used  for  determining  steady  state  response,  does  not  return  the  x  state 

PDSW.M-  MIMO  state  feedback  control  submarine  control  model,  with  wave  forces  and  trim 
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PDSWFF.M-  Same  as  PDSW.M,  with  disturbance  feedforward 

PDT.M-  used  for  determining  steady  state  response,  does  not  return  the  x  state 

Wave  Force  Data  files: 

For  the  first  order  motions  /  second  order  forces,  each  displacement  /  force  is  given  in 
terms  of  a  complex  number. 

Wave  spectral  data  files 

[wave  frequency  (rad/ sec),S(CO)  (),  amplitude  (feet),  phase  (rad)] 

SPEC_A.'FKT 

SPEC^-TXT 

SPEC_C.TXT 

SPEC_D.TKT 

First  order  motions 

[CO  (rad/sec),  COtencounter  (rad/sec),  surge,  sway,  heavy,  pitch,  yaw,  roU] 

LMOT_A.TXT 

LMOT_B.TKT 

LMOT_C.TXT 

LMOT_D.TXT 

Second  order  forces 

[COi  (rad/ sec),  (O2  (rad/ sec),  Fx,  Fy,  Fz,  Mx,  My,  Mz  ] 


FORC_A.TXT 

F0RC_B-'FXT 

FORC_C.TXT 

FORC_D.TXT 


118 


SUBOFF.M 


function  [A,B,MMI,Ettiat,  D,WS]— suboff(U); 

%  SDV  hydrodynamic  data,  TM  231-78  page  6 

ff=2;  %fudge  factor  to  make  planes  authority  realistic  for  300  ft  sub 

g=32.2;  % 

L=300;  %  feet 

p=1.94;  %  slug/ft^3; 


m  =  0.018296*0.5*p*L^3; 
Iz  =  0.001  084*0.5VL^5; 
ly  =  0.001 08*0.5  ; 

mxg  =-0.127467E-4; 

W=m*g; 

zgb=l;  %  feet 

xgb=0;  %  feet 


Mqdot  =-0.000860*0.5*p*L^5; 
Mwdot  =-0.000561*0.5*p*L^4; 

Mq  =-0.003702*0.5*p*L^4*U; 

Mw  =  0.010324*0.5*p*L^3*U; 

Mds  =-fP0.002409*0.5*p*L"'3*U"2; 
Mdb  =-Mds/4; 

Zqdot  =-0.000633*0.5VL^4; 

Zwdot  =-0.014529*0.5*p*L^3; 

Zq  =-0.007545*0.5*p*L^3*U; 

Zw  =-0.01 391 0*0.5*p*L'^2*U; 

Zds  =-fP0.005603*0.5*p*L''2*U^2; 
Zdb  =Zds/2; 


%  define  mass  matrix,  compute  mass  matrix  inverse 
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mm=eye(4); 
mm(l  ,1 ) =m-Zwdot; 
mm(l  ,2)=-Zqdot; 
mm  (2, 1 ) =-Mwdot; 
mm  (2,2) =Iy-Mqdot; 

A==zeros(4,4); 

A(l,l)=Zw; 

A(1 ,2)=m*U+Zq; 

A(2,l)=Mw; 

A(2,2)=Mq; 

A(2,3)=-2gb*W; 

A(3,2)=l; 

A(4,l)=l; 

A(4,3)=-U; 


B=[  Zdb  Zds;Mdb  Mds;0  0;0  0]; 


MMI=inv(mm); 


A=MMI*A; 

Emat=MMI(l  :2,1 :2)*[m*zgb,0;0,-m*zgb]; 
%  diameter  and  surface  area  calculation 
D=L/8.575; 


WS=67.651*(L/13.9792)"'2; 
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AXNL.M 


function  ax=axnl(x) 
global  Amat  Emat 
n=max(size(x)); 
a=Amat; 

a(4,l)=Amat(4,l)*cos(x(3)); 

x(3)=sin(x(3)); 

ax=a*x; 

ax(l  :2)=ax(l  :2)+Emat*[x(2)'^2pc(2)*x(l)]; 
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SMLQR.M 

function  [k,s,ks,e]=snilqr(a,b,q) 

%  [k,s,ks,e]=smlqr(a,b,c])  detemaines  the  sliding 
%  mode  control  law  for  the  system 

%  xdot=a*x+b*u  where 

%  u=-kx-kn*satsgn(sigma)  and 

%  sigma=s*x.  q  is  a  positive  definite 

%  symmetric  weighting  matrix 

%  used  in  assigning  error  weights 

%  (LQR)  to  the  states.  Uses  Utkin's  method 

%  as  detailed  in  Hawkinson  pplO-17 


[n,m]=size(b); 

%  do  transformation  if  required 
if  norm(b(m+l:n,l:ni))>eps‘^0.5 

[t,bl]=qr(b); 

bl=bl(l:m,l:m); 

t=t'; 

else 

t=eye(n); 

bl=b(l:m,l:m); 


q=t*q*t’; 

a=t*a*t'; 


qll=q(l:m,l:m); 
ql  2=q(l  :m,m+ 1  :n); 
q21=q(m+l:n,l:tn); 
q22=q(m+l  ;n,m+l  :n); 

all=a(l:m,l:m); 

al2=a(l:m,m+l:n); 

a21=a(m+l:n,l:tn); 

a22=a(m+l:n,m+l:n); 

as=a22-a21  *inv(ql  l)*ql  2; 
qs=q22-q21*inv(ql  l)*ql  2; 


kt=are(as,a21*inv(ql  l)*(a21')  ,qs); 
c2=(inv(ql  I)*(ql2+a21')*kt)'; 
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k=-rinv(bl)*(al  1  +c2'*a21),inv(bl)*(al2+c2'*a22)]* 
s=rref([eye(m);  c2]'*t)'; 


ks=-inv(s'*b); 


e=eig([a+  [bl  ;zeros(n-m,m)]*k*t']); 
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APPENDIX  B 


MAPLE®  Solutions 

Determination  of  MIMO  state  feedback  control  steady  state 
Determination  of  MIMO  sliding  mode  control  steady  state 


This  maple  script  determines  the  steady  state  response  in  the  vertical  plane 
using  MIMO  state  feedback  control.  Constant  disturbances  are  assumed  and  the 
linear  equations  of  motion  applied  to  find  the  steady  state  depth  error .  It 

is  also  assumed  that  both  control  deflections  are  less  than  maximum. _ 

EQl  and  EQ2  are  the  heave  and  pitch  linear  equations  of  motion _ 

>  EQl :  =  ( all*w*u+al3  *  theta+bll*u^2  *dl+bl2  *u^2*d2 )  +Fd; 

EQl  :=  all  w  u-\-  al3  Q  bl  1  dl  bl2  iP'  d2  +  Fd 

>  eQ2  :  =  (a21*w*u+a23*theta+b21*u''2*dl+b22*u^2*d2)  +Md; 

EQ2  :=  a21  w  u-\-  a23  0  +  b21  iP  dl  +  b22  iP  d2  -f  Md 

>  readlib ( isolate) : 


>  w  :=  theta*u; 

w  :=  0  M 

The  linear  state  feedback  laws  at  steady  state  are: 


>  dl  :=  Kll*w+K13*theta+K14*z; 

dl  :=K11  Qu  +  K13Q  +  K14z _ 

>  d2  :=  K21*w+K23*theta+K14*z; 

d2-K21Qu  +  K23d-^K14z _ 

After  the  application  of  the  control  laws,  the  equations  of  heavy  and  pitch 
become : 

>  EQ1=0; 

all  Qu‘^  +  al3e  +  bll  {Kll  Qu  +  K13%  +  K14  z)  +  bl2  u^{K21du  +  K23  d  +  KMz) 

+  Fd  =  0 _ _ _ _ _ 

>  EQ2=0; 

a21Qu^  +  a23Q  +  b21  {K1 1  Q  u  + K13  Q  +  K14  z)  +  b22  {K21  %  u  + K23  Q  +  K14  z) 

+  Md  =  Q _ _ _ _ 

Remove  z  from  EQl  and  EQ2  _ _ _ _ _ — 

>  FI :  =coef f (collect {expand(EQl)  , z)  , z, 1) ; 

FI  :=bl2u^  K14  +  bll  KM  _ _ 

>  F2 : =coef f (collect ( expand (EQ2 ) , z) , z , 1) ; 

F2  :=  b22  KM  +  b21  KM  _ 

>  EQ3:=siinplify(EQl-EQ2*Fl/F2)  ; 

EQ3  :=  (a77  9  b22  +  all  0  b21  +  bll  Kll  0  b22  +  bll  K13  9  b22  +  al3  0  b22 

+  al3Qb21  -a23Bbll  +Fd  b22  +  Fd  b21  -  Md  bl2  -Mdbll+  bl2  K21  9  b21 
+  bl2  K23  Qb21  -a21  d  bl2  -  a21  0  bll  -  b21  Kll  0  bl2  -  b21  K13  0  bl2 

-b22u^  K21  Bbll-  b22  K23  Bbll  -a23Q  bl2)l {b22  +  b21) _ _ 

>  isolate (EQ3 , theta) ; 

0  =  {-Fdb22  -  Fdb21  +  Md  bl2  +  Md  bl})l{all  b22  +  all  b21  +  bll  Kll  b22 
+  bll  u'^  K13b22  +  al3  b22  +  al3  b21  -  a23  bll  +  bl2  K21  b21  +  bl2  K23  b21 
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-  all  u-  bl2  -  all  bll  -  bll  Kll  bll  -  bll  u'^  K13  bll  -  bll  Kll  bll 
-bllu^KlS  bll-al3bll) 

Find  the  steady  state  value  of  theta  by  setting  EQ3=Q 

>  thetas s  :  =:coef  f  (collect  (EQ3  ,  theta)  ,  theta,  0)  /coef  f  (collect  (  EQ3  ,  theta)  ,  theta,  1 )  ; 

thetass  ;=  {Fdbll  +  Fdbll  -  Md  bll  -  Md  bll)l{all  bll  +  all  bll 

+  bll  Kll  bll  +  bll  K13  bll  +  al3  bll  +  al3  bll  -  al3  bll  +  bll  Kll  bll 

+  bll  u'^  K13  bll  -  all  bll  -  all  \?  bll  -  bll  Kll  bll  -  bll  K13  bll 

-  bll  Kll  bll  -  bll  K13  bll  -  al3  bll) 

Substitute  thetass  into  EQl  and  solve  for  steady  state  z _ 

>  temp: =coef f (collect (EQl, z) , z , 0) /coef f (collect (EQl , z) , z , 1 ) ; 

all  9  e  +  bll  {KllBu  +  K13Q)  +  bll  {Kll  Q  u  + K13  d)  +  Fd 

bllu^  K14  + bll  u^KM 

>  zss 

>  coef  f  (collect  (temp,  theta)  ,  theta,  1)  *  thetas  s  + coef  f  (collect  { temp,  theta)  ,  theta,  0) 

>  ; 

2^5  ■.=  {al3  +  bll  {Kll  u  +  K13)  +  bll  {Kll  u  +  K13)  +  al  1  u^) 

{Fd  bll  +  Fd  bll  -  Md  bll  -  Md  bll)  l((bl2  K14  +  bll  Km)  {all  bll 

+  all  bll  +  bll  Kll  bll  +  bll  i?  K13  bll  +  al3  bll  +  al3  bll  -  al3  bll 

+  bll  Kll  bll  +  bll  K13  bll  -  all  bll  -  all  u^bll  -  bll  Kll  bll 

-  bll  K13  bll  -  bll  Kll  bll  -  bll  K13  bll  -  al3  bll)) 

Fd 

_| - 

bllu^  K14  +  bll  u'^  K14 

Determine  the  coefficients  Cl  and  C2  such  that  zss=Cl *Fd+C2 *Md 

>  Cl : ^simplify (coef f (collect (expand (zss) , Fd) , Fd, 1 ) ) ; 

Cl  -.^{lall  bll  +  2  all  bll  -  all  bll  -  all  u^bll  -  al3  bll  +  2  al3  bll 

-  al3  bll  +  2  al3  bll  +  2  bll  Kll  bll  +  2  bll  K13  bll  -  bll  Kll  bll 
+  2  bll  Kll  bll  +  2  bll  K13  bll  -  bll  Kll  bll  -  bll  K13  bll 

-  bll  K13  bll  +  bll  K13  bll  +  bll  Kll  bll  +  bll  Kll  bll  +  bll  K13  bll) 

/(m^  K14  {bll  +  bll )  {all  bll  +  all  bll  +  bll  Kll  bll  +  bll  u'^  K13  bll 

+  al3  bll  +  al3  bll  -  al3  bll  +  bll  Kll  bll  +  bll  K13  bll  -  all  bll 

-  all  bll  -  bll  V?  Kll  bll  -  bll  i?  K13  bll  -  bll  Kll  bll  -  bll  K13  bll 
-al3  bll)) 

>  C2 : =simpli fy ( coef f ( collect (expand (zss) , Md) , Md, 1 ) ) ; 

Cl  :=  -  {al3  +  bll  Kll  +  bll  K13  +  bll  Kll  +  bll  K13  +  all  w^)/(( 
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> 


al  1  b22  +  all  b21  +  bll  Kll  b22  +  bll  i?  K13  b22  +  al3  b22  +  al3  b21  -  a23  bll 

+  bl2u^  K21  b21  +  bl2  K23  b21  -  a21  bl2  -  a21  bll  -  b21  Kll  bl2 

-  b21  K13  bl2  -  b22  K21  bll  -  b22  K23  bll  -  a23  bn)  Km) _ 

Check  that  zss=Cl*Fd+C2*Md _ _ _ _ _ 

eq6 : =simplify(zss-Cl*Fd-C2*Md)  ; 

eq6  0 
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This  maple  script  determines  the  steady  state  response  in  the  vertical  plane 
with  basic  MIMO  sliding  mode  control.  Constant  disturbances  are  assumed  and 
the  linear  equations  of  motion  applied  to  find  the  steady  state  depth  error. 
It  is  also  assumed  that  both  control  deflections  are  less  than  maximum  and 

that  neither  sliding  surface  is  saturated. _ 

EQl  and  EQ2  are  the  heave  and  pitch  linear  equations  of  motion _ 

>  EQl :  =  (all*w*u+al3*  theta+bll*u'^2*dl+bl2*u'^2*d2  )  +Fd; 

EQl  ;=  allwu-^ al3  0  +  bll  dl  +  bl2  \?  d2  +  Fd 

>  eQ2  ; = {a21*w*u+a23*theta+b21*u^2*dl+b22*u^2*d2 ) +Md; 

EQ2  :=  a2]  w  w  +  a23  0  +  b21  dl  +  b22  iP'  d2  +  Md 

>  readlib { isolate)  : 


>  w  :=  theta*u; 


w  :=  0  w 


The  linear  state  feedback  laws  at  steady  state  are: 


>  dl  :=  Kll*w+K13*theta+Ksll*etal*sigmal/phil+Ksl2*eta2*sigma2/phi2; 

phil  phi2 


>  d2  :=  K21*w+K23*theta+Ks21*etal*siginal/phil+Ks22*eta2*sigma2/phi2; 

^  ^  Ks21  etal  sigmal  Ks22  eta2  sigma2 

d2  :=K21  Qu  +  K23%-\- - -rr^ - + - - 

phil  phi2 

>  sigitial;=Sll*w+S13*theta+S14*z; 

sigmal  :=  Sll  Q  u  +  S13  0  +  S14  z 


>  sigina2  :=S21*w+S23*theta+S24*2; 

j■^gma2  :=  527  0  M  +  S23  0  +  S24  z 


After  the  application  of  the  control  laws,  the  equations  of  heavy  and  pitch 
become : 


>  EQ1=0; 


all  0  +  al3Q  +  bll 


( 


/i:77  0  M  +  7i:7i  0 + 


V 


Ksll  etal  (577  0  M +  573  0  + 574z) 
phil 


+  ■ 


+  bl2u^  K21BU  +  K23Q 


Ksl2  eta2  (527  0  h  +  52J  0  +  S24  z) 
phi2 

Ks21  etal  (577  0  «  +  575  0  +  S14  z)  ^  Ks22  eta2  (S21  0  h  +  525  0  +  524  z)  ^  ^  ^ 

phil  phi2 


>  EQ2=0; 

9  .  Ksll  etal  (Sll  Qu  +  S13Q  +  S14z) 

a21  Qu^  +  a23Q  +  b21  u^\K11Qu  +  K13  0  + - — - 


+ 


Ksl2  eta2  (527  0  m  +  525  0  +  524  z) 
phi2 


+  b22 


K21Bu  +  K23  0 


Ks21  etal  (Sll  B  u  +  575  0  +  57^  z)  ^  Ks22  eta2  {S21  0  m  +  523  0  +  S24  z )  ^  ^  _  q 

phil  ^  phi2  )  _ 


Remove  z  from  EQl  and  EQ2 _ _ _ 

>  FI ; =coef f (collect (expand (EQl) , z ) , z , 1 ) ; 
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bll  Ksl2  eta2  S24  bll  Ksll  etal  S14  ^  bl2  Ks22  eta2  S24 

phil  phi2 

bll  Ks21  etal  S14 
^  phil 

>  F2 : =coeff (collect ( expand (EQ2)  ,  z)  ,  z,  1)  ; 

bll  Ksll  etal  S14  bll  Ksll  etal  S14  bll  Ksll  etal  S24 

FI  := - + - tt; - + - 7^ - 

phi2  phil  phi2 

bll  Ksll  etal  S14 

phil  _ 

>  EQ3:=simplify{EQl-EQ2*Fl/F2)  ; 

EQ3  :=  -  {-bll  Ksll  etal  Sll  6  bll  Ksll  etal  S14  -  bll  Kll  d  phil  bll  Ksll  etal  S24 
-bll  Ksll  etal  S13  0  bll  Ksll  etal  S24  -bll  Ksll  etal  Sll  0  bll  Ksll  etal  S14 

-bll  Ksll  etal  S23  0  bll  Ksll  etal  S14  -bll  Kll  0  phil  bll  Ksll  etal  S14 

-bll  K13  0  phil  bll  Ksll  etal  S14  -bll  K13  0  phil  bll  Ksll  etal  S14 
-all  Q  u'^  phil  bll  Ksll  etal  S24-all  0  phil  bll  Ksll  etal  S14 

-all  phil  bll  Ksll  etal  S14  -allQu^  phil  bll  Ksll  etal  S14 

-  al3  0  phil  bll  Ksll  etal  S14  -  al3  0  phil  bll  Ksll  etal  S14 

-  al3  0  phil  bll  Ksll  etal  S14  -  al3Q  phil  bll  Ksll  etal  SI 4 

-bllu^  Kll  Q  phil  bll  Ksll  etal  S24  -  bll  Kll  0  phil  bll  Ksll  etal  S14 

-  bll  K13  0  phil  bll  Ksll  etal  S24  -  bll  K13  0  phil  bll  Ksl  1  etal  SI 4 

-bllu^  Ksll  etal  Sll  6  bll  Ksll  etal  S14  -  bll  Ksll  etal  S13  0  bll  Ksll  etal  S14 

-bllu^  Ksll  etal  Sll  0  bll  Ksll  etal  S14  -  bll  Ksll  etal  S13  0  bll  Ksll  etal  S14 

-  Fd phil  bll  Ksll  etal  S24  -  Fd phil  bll  Ksll  etal  S14  -  Fd phil  bll  Ksll  etal  S14 

-  Fd  phil  bll  Ksll  etal  S14  +  all  0  phil  bll  Ksll  etal  S24 

+  a21  Q  phil  bll  Ksll  etal  S14  +  all  0  u'^  phil  bll  Ksll  etal  S24 

+  all  0  phil  bll  Ksll  etal  S14  +  al3  0 phil  bll  Ksll  etal  S24 
+  a23  0  phil  bll  Ksll  etal  S14  +  al3  0  phil  bll  Ksll  etal  S24 

+  a23  0  phil  bll  Ksll  etal  S14  +  bll  Kll  0  phil  bll  Ksll  etal  S24 

+  bll  Kll  Q  phil  bll  Ksll  etal  S14  +  bll  K13  0  phil  bll  Ksll  etal  S14 

+  bll  K13  0  phil  bll  Ksll  etal  S14  +  bll  Ksll  etal  Sll  Q  bll  Ksll  etal  S14 

+  bll  Ksl  1  etal  SI 3  0  bll  Ksll  etal  S24  +  bll  Ksll  etal  Sll  0  bll  Ksll  etal  S14 

+  bll  Ksll  etal  S13  0  bll  Ksll  etal  S14  +  bll  Kll  B  phil  bll  Ksll  etal  S14 

+  bll  Kll  B  phil  bll  Ksll  etal  S14  +  bl2  K23  B  phil  bll  Ksll  etal  S14 
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+  b22  K23  e  phi2  bll  Ksll  etal  SI4  +  b22  Ks21  etal  SllQbll  Ksl2  eta2  S24 

+  b22  Ks21  etal  SI 3  0  bll  Ksl2  eta2  S24  +  b22  Ks22  eta2  S21  0  bll  Ksll  etal  SI 4 

+  b22  Ks22  eta2  S23  0  bll  Ksll  etal  S14  +  Mdphil  bll  Ksl2  eta2  S24 

+  Md phi2  bll  Ksll  etal  S14  +  Md phil  bl2  Ks22  eta2  S24  +  Md phi2  bl2  Ks21  etal  S14)/( 
b21  Ksl2  eta2  S24  phil  +  b21  Ksll  etal  S14  phi2  +  b22  Ks22  eta2  S24  phil 
+  b22Ks21  etal  S14 phi2) 

>  isolate {EQ3 , theta) ; 

B  =  {Fd phil  b21  Ksl2  eta2  S24  +  Fd phi2  b21  Ksll  etal  S14  +  Fd phil  b22  Ks22  eta2  S24 
+  Fd  phi2  b22  Ks21  etal  S14  -  Md  phil  bll  Ksl2  eta2  S24  -  Md  phi2  bll  Ksll  etal  S14 

-  Md  phil  bl2  Ks22  eta2  S24  -  Md  phi2  bl2  Ks21  etal  S14)l( 

-bl  1  Ksll  etal  Sll  b22  Ks22  eta2  S24  -bll  Kll  phil  b22  Ks22  eta2  S24 
-bll  Ksl2  eta2  S21  b22  Ks21  etal  S14  -bll  Ksl2  eta2  S23  b22  Ks21  etal  S14 
-all  iP'  phil  b21  Ksl2  eta2  S24  -  bll  P  Ksll  etal  SI  3  b22  Ks22  eta2  S24 
-bll  P  Kll  phi2  b22  Ks21  etal  S14  -  bll  P  K13  phil  b22  Ks22  eta2  S24 

-  bll  P  K13  phi2  b22  Ks21  etal  S14  -  al3  phi2  b21  Ksll  etal  S14 

-  all  P  phi2  b21  Ksll  etal  S14  -all  P  phil  b22  Ks22  eta2  S24 
-all  P  phi2  b22  Ks21  etal  S14  -  al3  phil  b21  Ksl2  eta2  S24 
-bl2p  K23  phi2  b21  Ksll  etal  S14  -  al3  phil  b22  Ks22  eta2  S24 

-  al3  phi2  b22  Ks21  etal  S14  -  bl2  P  K21  phil  b21  Ksl2  eta2  S24 
-bl2p  K21  phi2  b21  Ksll  etal  S14  -  bl2  P  K23 phil  b21  Ksl2  eta2  S24 
+  a21  P  phi2  bl2  Ks21  etal  S14  +  a23  phil  bll  Ksl2  eta2  S24 

+  a23  phi2  bll  Ksll  etal  S14  -  bl2  P  Ks21  etal  Sll  b21  Ksl2  eta2  S24 

-  b  12  P  Ks21  etal  SI 3  b21  Ksl2  eta2  S24  -  bl2  P  Ks22  eta2  S21  b21  Ksll  etal  S14 

-  bl2  P  Ks22  eta2  S23  b21  Ksll  etal  S14  +  a21  P  phil  bll  Ksl2  eta2  S24 
+  a21  P  phi2  bll  Ksll  etal  S14  +  a21  P  phil  bl2  Ks22  eta2  S24 

+  b21  P  Ksl2  eta2  S21  bl2  Ks21  etal  S14  +  b21  P  Ksl2  eta2  S23  bl2  Ks21  etal  S14 

+  b22  P  K21  phil  bll  Ksl2  eta2  S24  +  b22  P  K21  phi2  bll  Ksll  etal  S14 
+  a23  phil  bl2  Ks22  eta2  S24  +  a23  phi2  bl2  Ks21  etal  S14 

+  b21  P  Kll  phil  bl2  Ks22  eta2  S24  +  b21  P  Kll  phi2  bl2  Ks21  etal  S14 

+  b21  P  K13  phil  bl2  Ks22  eta2  S24  +  b21  P  K13  phi2  bl2  Ks21  etal  S14 

+  b21  P  Ksll  etal  Sll  bl2  Ks22  eta2  S24  +  b21  P  Ksll  etal  S13  bl2  Ks22  eta2  S24 

+  b22  P  K23  phil  bll  Ksl2  eta2  S24  +  b22  P  K23  phi2  bll  Ksll  etal  S14 
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+  b22  Ks21  etal  Sll  bll  Ksl2  eta2  S24  +  b22  Ks21  etal  S13  bll  Ksl2  eta2  S24 
+  b22  Ks22  eta2  S21  bll  Ksll  etal  S14  +  b22  Ks22  eta2  S23  bll  Ksll  etal  Sm) 

Find  the  steady  state  value  of  theta  by  setting  EQ3  =  0 _ 

>  thetass  :  =siinplify  (coef  f  (collect  (EQ3  ,  theta)  ,  theta,  0 )  /coef  f  (collect  (EQ3  ,  theta)  , 

>  theta, 1) ) ; 

thetass  :=  {Md phi2  bll  Ksll  etal  S14  +  Mdphil  bl2  Ks22  eta2  S24 
+  Md phi2  bl2  Ks21  etal  S14  -  Fd phil  b21  Ksl2  eta2  S24  +  Md phil  bll  Ksl2  eta2  S24 

-  Fd  phil  b22  Ks22  eta2  S24  -  Fd  phi2  b22  Ks21  etal  S14  -  Fd  phi2  b21  Ksll  etal  S14)l{ 
-bll  Ksll  etal  Sll  b22  Ks22  eta2  S24  -  bll  Kll  phil  b22  Ks22  eta2  S24 

-bll  Ksl2  eta2  S21  b22  Ks21  etal  S14  -  bll  Ksl2  eta2  S23  b22  Ks21  etal  S14 

-  all  phil  b21  Ksl2  eta2  S24  -  bll  Ksll  etal  S13  b22  Ks22  eta2  S24 

-bll  Kll  phi2  b22  Ks21  etal  S14  -bll  K13  phil  b22  Ks22  eta2  S24 

-  bll  K13  phi2  b22  Ks21  etal  S14  -  al3  phi2  b21  Ksll  etal  S14 

-  all  iP"  phi2  b21  Ksll  etal  S14  —  all  iP  phil  b22  Ks22  eta2  S24 
-all  u'^  phi2  b22  Ks21  etal  S14  -  al3  phil  b21  Ksl2  eta2  S24 
-bl2u^  K23  phi2  b21  Ksll  etal  S14  -  al3  phil  b22  Ks22  eta2  S24 

-  al3  phi2  b22  Ks21  etal  S14  -  bl2  K21  phil  b21  Ksl2  eta2  S24 

-  bl2  K21  phi2  b21  Ksll  etal  S14  -  bl2  iP  K23  phil  b21  Ksl2  eta2  S24 

+  a21  phi2  bl2  Ks21  etal  S14  +  a23  phil  bll  Ksl2  eta2  S24 

+  0.23  phi2  bll  Ksll  etal  S14  -  bl2  Ks21  etal  Sll  b21  Ksl2  eta2  S24 
-bl2u^  Ks21  etal  SI 3  b21  Ksl2  eta2  S24  -  bl2  Ks22  eta2  S21  b21  Ksll  etal  S14 

-  bl2  iP  Ks22  eta2  S23  b21  Ksll  etal  S14  +  a21  iP phil  bll  Ksl2  eta2  S24 

+  a21  phi2  bll  Ksll  etal  S14  +  a21  phil  bl2  Ks22  eta2  S24 

+  b21  Ksl2  eta2  S21  bl2  Ks21  etal  S14  +  b21  Ksl2  eta2  S23  bl2  Ks21  etal  S14 

+  b22  K21  phil  bll  Ksl2  eta2  S24  +  b22  K21  phi2  bll  Ksll  etal  S14 

+  a23  phil  bl2  Ks22  eta2  S24  +  a23  phi2  bl2  Ks21  etal  S14 

+  b21  Kll  phil  bl2  Ks22  eta2  S24  +  b21  P’  Kll  phi2  bl2  Ks21  etal  S14 

+  b21  P  K13  phil  bl2  Ks22  eta2  S24  +  b21  P  K13  phi2  bl2  Ks21  etal  S14 

+  b21  P  Ksll  etal  Sll  bl2  Ks22  eta2  S24  +  b21  P  Ksll  etal  S13  bl2  Ks22  eta2  S24 

+  b22  P  K23  phil  bll  Ksl2  eta2  S24  +  b22  P  K23  phi2  bll  Ksll  etal  S14 

+  b22  P  Ks21  etal  Sll  bll  Ksl2  eta2  S24  +  b22  P  Ks21  etal  S13  bll  Ksl2  eta2  S24 

+  b22  P  Ks22  eta2  S21  bll  Ksll  etal  S14  +  b22  P  Ks22  eta2  S23  bll  Ksll  etal  Sm') 

^Substitute  thetass  into  EQl  and  solve  for  steady  state  z _ _ _ 
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>  temp : =coef f (collect  (EQl , z) , z , 0 ) /coef f (collect (EQl ,  z)  , z , 1 )  ; 
f  o 


temp  := 


all  Qu^  +  aI3Q 


21  Ksll  etal  (Sll  du  +  SlSe)  Ksl2  eta2  (S21  6  u  +  S23  6) 

+  bll  Kll  e  u  +  K13e+ - - ^  +  -  ^ 


phil 


phi2 


21  ^  Ks21etal(Slldu  +  S13e)  Ks22  eta2  (S21  6  u  +  S23  6) 

+  bl2  K21  eu  +  K23e+ - r ^  + 


phil 


+  Fd 


bll 


V 

Ksl  1  etal  S14  Ksl2  eta2  S24 


phil 


phi2 


+  bl2 


J 


phi2 


Ks21  etal  S14  Ks22  eta2  S24 


phil 


■  +  ■ 


phi2 


>  zss  := 

>  coeff (collect (temp, theta) , theta, 1) *thetass+coef f (collect (temp,  theta) , theta, 0) 


f 


ZSS  := 


al3  +  bll 


fr..,  .  r.,.  .  Ksll  etal  (Sll  U  +  S13)  .  Ksl2  eta2  (S21  u  +  S23)] 

K1 1  u  +  K13  ^ - T-r,  I"  , 

phil  phi2 


J 

.b,2udK2,u.  K23  .  1 .  uA  ( 

phil  phi2  J  J 

Md phi2  bll  Ksll  etal  S14  +  Md phil  bl2  Ks22  eta2  S24  +  Md phi2  bl2  Ks21  etal  S14 

-  Fd  phil  b21  Ksl2  eta2  S24  +  Md  phil  bll  Ksl2  eta2  S24  -  Fd  phil  b22  Ks22  eta2  S24 


Fd  phi2  b22  Ks21  etal  S14  -  Fd  phil  b21  Ksll  etal  S14) 


r 


bll  u 


V 


Ksl  1  etal  S14  Ksl2  eta2  S24 


phil 


phi2 


+  bl2u' 


Ks21  etal  S14  Ks22  eta2  S24 
-  +  - 


phil 


phi2 


-bll  Ksll  etal  Sll  b22  Ks22  eta2  S24  -  bll  Kll  phil  b22  Ks22  eta2  S24 
-bll  Ksl  2  eta2  S21  b22  Ks21  etal  S14  -bll  Ksl  2  eta2  S23  b22  Ks21  etal  S14 

-all  phil  b21  Ksl2  eta2  S24  -  bll  Ksll  etal  S13  b22  Ks22  eta2  S24 

-bll  Kll  phi2  b22  Ks21  etal  S14  -bll  K13  phil  b22  Ks22  eta2  S24 

-bll  K13  phi2  b22  Ks21  etal  S14  -  al3  phi2  b21  Ksll  etal  S14 

-all  phi2  b21  Ksll  etal  S14  -all  phil  b22  Ks22  eta2  S24 

-all  phi2  b22  Ks21  etal  S14  -  al3  phil  b21  Ksl  2  eta2  S24 

-bl2u^  K23  phi2  b21  Ksll  etal  S14  -  al3  phil  b22  Ks22  eta2  S24 

-  al3  phi2  b22  Ks21  etal  S14  -  bl2  K21  phil  b21  Ksl2  eta2  S24 

-  bl2  K21  phi2  b21  Ksll  etal  S14  -  bl2  K23  phil  b21  Ksl2  eta2  S24 

+  a21  phi2  bl2  Ks21  etal  S14  +  a23  phil  bll  Ksl2  eta2  S24 

+  a23  phi2  bll  Ksll  etal  S14  -  bl2  Ks21  etal  Sll  b21  Ksl2  eta2  S24 

-  bl2  Ks21  etal  S13  b21  Ksl2  eta2  S24-  bl2  Ks22  eta2  S21  b21  Ksll  etal  S14 
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-  bl2  u-  Ks22  eta2  S23  b21  Ksll  etal  S14  +  a2]  phi  1  bll  Ksl2  eta2  S24 

+  a21  phi2  bll  Ksll  etal  SI 4  +  a21  phil  bl2  Ks22  eta2  S24 

+  b21  Ksl2  eta2  S21  bl2  Ks21  etal  S14  +  b21  Ksl2  eta2  S23  bl2  Ks21  etal  S14 

+  b22  K21  phil  bll  Ksl2  eta2  S24  +  b22  K21  phi2  bll  Ksll  etal  S14 

+  a23  phil  bl2  Ks22  eta2  S24  +  a23  phi2  bl2  Ks21  etal  S14 

+  b21  Kll  phil  bl2  Ks22  eta2  S24  +  b21  Kll  phi2  bl2  Ks21  etal  S14 

+  b21  K13  phil  bl2  Ks22  eta2  S24  +  b21  K13  phi2  bl2  Ks21  etal  S14 

+  b21  Ksll  etal  Sll  bl2  Ks22  eta2  S24  +  b21  Ksll  etal  S13  bl2  Ks22  eta2  S24 

+  b22  K23  phil  bll  Ksl2  eta2  S24  +  b22  K23 phi2  bll  Ksll  etal  S14 

+  b22  Ks21  etal  Sll  bll  Ksl2  eta2  S24  +  b22  Ks21  etal  S13  bll  Ksl2  eta2  S24 

+  b22  Ks22  eta2  S21  bll  Ksll  etal  S14  +  b22  Ks22  eta2  S23  bll  Ksll  etal  Sm) 


> 


Fd 


+  • 


bll  u 


Ksll  etal  S14  Ksl2  eta2  S24  '\ 


phil 


phi2 


j  +  bl2u^ 


Ks21  etal  S14  Ks22  etal  S24 


phil 


phi2 


Determine  the  coefficients  Cl  and  C2  such  that  zss-Cl*Fd+C2 *Md 


Cl :  =  simplify (coeff (collect (expand (zss)  , Fd) , Fd, 1) )  ; 


Cl  -  (2  phl'^phil  al3  b22  Ks21  etal  S14  +  2  phi2  phil^  al3  b21  Ksl2  eta2  S24 
H-  2  phi2  phil^  al3  b22  Ks22  eta2  S24  +  2  phi2^  phil  al3  b21  Ksll  etal  S14 
+  phi2  phil^  bll  Kll  b21  Ksl2  eta2  S24  +  2  phi2  phil^  bl  1  u  Kll  b22  Ks22  eta2  S24 

+  phi2  phil  bll  Ksll  etal  S13  b21  Ksl2  eta2  S24 

+  2  phi2  phil  bll  Ksll  etal  SI 3  b22  Ks22  eta2  S24 

+  2  phi2  phil  bll  Ksll  etal  Sll  b22  Ks22  eta2  S24 

+  2  phi2^  phil  bll  K13  b22  Ks21  etal  S14  +  phi2^ phil  bll  K13  b21  Ksll  etal  S14 

+  2  phi2  phil^  bll  K13  b22  Ks22  eta2  S24  +  phi2'^ phil  bll  Kll  b21  Ksll  etal  S14 

H- phi2  phil  bl  1  Ksll  etal  Sll  b21  Ksl2  eta2  S24 

+  phi2  phil^  bll  K13  b21  Ksl2  eta2  S24  +  2  phi2^ phil  bll  Kll  b22  Ks21  etal  S14 

+  2  phi2  phil  bll  Ksl2  eta2  S21  b22  Ks21  etal  S14 

+  phi2  phil  bll  Ksl2  eta2  S21  b21  Ksll  etal  S14 

+  phi2  phil  bll  Ksl2  eta2  S23  b21  Ksll  etal  S14 

+  2  phi2  phil  bll  Ksl2  eta2  S23  b22  Ks21  etal  S14 
+  2  phi2  phil^  bl2  K21  b21  Ksl2  eta2  S24 
+  2  phi2  phil  bl2  Ks21  etal  Sll  b21  Ksl2  eta2  S24 
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+  phi2  phil^  bl2  K23  b22  Ks22  eta2  S24  +  2  phi2  phil^  bl2  K23  b21  Ksl2  eta2  S24 

+  2  phi2^phil  bl2  K21  b21  Ksll  etal  S14  +  pha'^  phil  bl2  K21  b22  Ks21  etal  S14 

+  phi2  phil^  bl2  K21  b22  Ks22  eta2  S24  +  2  pha'^  phil  bl2  K23  b21  Ksll  etal  S14 

+  phi2^  phil  bl2  K23  b22  Ks21  etal  SI 4 
+  phi2  phil  bl2  Ks21  etal  S13  b22  Ks22  eta2  S24 
+  2  phi2  phil  bl2  Ks21  etal  S13  b21  Ksl2  eta2  S24 
+  phi2  phil  bl2  Ks21  etal  Sll  b22  Ks22  eta2  S24 
+  phi2  phil  bl2  Ks22  eta2  S21  b22  Ks21  etal  S14 
+  2  phi2  phil  bl2  Ks22  eta2  S21  b21  Ksll  etal  S14 
+  phi2  phil  bl2  Ks22  eta2  S23  b22  Ks21  etal  S14 
+  2  phi2  phil  bl2  iP"  Ks22  etal  S23  b21  Ksll  etal  SI 4 

+  2  phi2  phil^all  b21  Ksl2  eta2  S24  +  2 phi2  phil^  all  b22  Ks22  eta2  S24 
+  2  phi2^  phil  all  b22  Ks21  etal  S14  +  2 phi2^ phil  all  b21  Ksll  etal  S14 
+  phi2^  bll  Ksll  etal^  Sll  b22  Ks21  S14  +  phi2^  bll  Ksll^  etal^  Sll  b21  S14 

+  phil^  bll  Ksl2'^  eta2'^  S23  b21  S24  +  phil^bll  Ksl2  eta2'^  S23  b22  Ks22  S24 

+  pha'^  bl  1  Ksll  etal^  S13  b22  Ks21  S14  +  phi2^  bll  Ksll^  etal^  S13  b21  S14 

+  phil^  bll  Ksl2^  eta2'^  S21  b21  S24  +  phil'^  bll  Ksl2  eta2^  S21  b22  Ks22  S24 
+  phi2^  bl2  Ks21^  etal^  S13  b22  S14  +  p/z/2^  bl2  Ks21  etal^  S13  b21  Ksll  S14 
+  phil^  bl2  Ks22  eta2^  S21  b21  Ksl2  S24  +  phil^  bl2  i?  Ks22^  eta2^  S21  b22  S24 

+  phi2^  bl2  Ks2l'^  etal^  Sll  b22  S14  +  phi2^  bl2  Ks21  etal^  Sll  b21  Ksll  S14 

+  phil^  bl2  iP  Ks22  eta2^  S23  b21  Ksl2  S24  +  phil^  bl2  Ks22'^  eta2^  S23  b22  S24 

-  phi2^  phil  a21  bl2  Ks21  etal  S14-phi2  phil^  a23  bll  Ksl2  eta2  S24 

-  phi2^  phil  a23  bll  Ksll  etal  S14  -  phi2  phil^  a21  iP  bll  Ksl2  eta2  S24 

-  phi2^  phil  a21  P  bll  Ksll  etal  S14  - phi2  phiP  a21  P  bl2  Ks22  eta2  S24 

-  phi2  phil  b21  P  Ksl2  eta2  S21  bl2  Ks21  etal  S14 

-  phi2  phil  b21  P  Ksl2  eta2  S23  bl2  Ks21  etal  S14 

-phi2  phiP  b22  P  K21  bll  Ksl2  eta2  S24  - phi2^ phil  b22  P  K21  bll  Ksll  etal  S14 

-  phi2  phiP  a23  bl2  Ks22  eta2  S24  -  phi2^  phil  a23  bl2  Ks21  etal  S14 

-  phi2  phiP  b21  P  Kll  bl2  Ks22  eta2  S24  -  phi2^  phil  b21  P  Kll  bl2  Ks21  etal  S14 

-  phi2  phiP  b21  P  K13  bl2  Ks22  eta2  S24  -  phi2^  phil  b21  P  K13  bl2  Ks21  etal  S14 
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-phi2  phil  b21  Ksll  etal  Sll  bl2  Ks22  eta2  S24 

-  phi2  phil  b21  Ksll  etal  SI 3  bl2  Ks22  eta2  S24 

-  phi2  phil^  b22  K23  bll  Ksl2  eta2  S24  -  phi2^  phil  b22  K23  bl  1  Ksll  etal  S14 

-  phi2  phil  b22  Ks21  etal  Sll  bll  Ksl2  eta2  S24 

-  phi2  phil  b22  Ks21  etal  S13  bll  Ksl2  eta2  S24 

-  phi2  phil  b22  Ks22  eta2  S21  bll  Ksll  etal  S14 

-phi2  phil  b22  Ks22  eta2  S23  bll  Ksll  etal  SM^jiu'^  (phil  bll  Ksl2  eta2  S24 

+  phi2  bll  Ksll  etal  S14  +  phil  bl2  Ks22  eta2  S24  +  phi2  bl2  Ks21  etal  S14)  ( 

-bll  Ksll  etal  Sll  b22  Ks22  eta2  S24  -bll  Kll  phil  b22  Ks22  eta2  S24 
-bll  Ks  12  eta2  S21  b22  Ks21  etal  SI  4  -bll  Ksl2  eta2  S23  b22  Ks21  etal  SI  4 

-all  u'^  phil  b21  Ksl2  eta2  S24  -bll  Ksll  etal  SI  3  b22  Ks22  eta2  S24 
-bll  Kll  phi2  b22  Ks21  etal  S14  -  bll  K13  phil  b22  Ks22  eta2  S24 
-bll  K13  phi2  b22  Ks21  etal  S14  -  al3  phi2  b21  Ksll  etal  S14 
-all  phi2  b21  Ksll  etal  S14  -all  phil  b22  Ks22  eta2  S24 
-all  u'^  phi2  b22  Ks21  etal  S14  -  al3  phil  b21  Ksl2  eta2  S24 

-  bl2  K23  phi2  b21  Ksll  etal  S14  -  al3  phil  b22  Ks22  eta2  S24 

-  al3  phi2  b22  Ks21  etal  S14  -  bl2  K21  phil  b21  Ksl2  eta2  S24 

-  bl2  K21  phi2  b21  Ksll  etal  S14  -  bl2  K23  phil  b21  Ksl2  eta2  S24 
+  a21  iP'  phi2  bl2  Ks21  etal  S14  +  a23  phil  bll  Ksl2  eta2  S24 

+  a23  phi2  bll  Ksll  etal  S14  -  bl2  P  Ks21  etal  Sll  b21  Ksl2  eta2  S24 

-  bl2  P  Ks21  etal  S13  b21  Ksl2  eta2  S24  -  bl2  P  Ks22  eta2  S21  b21  Ksll  etal  S14 
-bl2  P  Ks22  eta2  S23  b21  Ksll  etal  S14  +  a21  P  phil  bll  Ksl2  eta2  S24 

+  a21  P  phi2  bll  Ksll  etal  SI 4  +  a21  p-  phil  bl2  Ks22  eta2  S24 

-\-b21  P  Ksl2  eta2  S21  bl2  Ks21  etal  S14  +  b21  P  Ksl2  eta2  S23  bl2  Ks21  etal  S14 

+  b22  P  K21  phil  bll  Ksl2  eta2  S24  +  b22  P  K21  phi2  bll  Ksll  etal  SI 4 
+  a23  phil  bll  Ksll  etal  S14  +  al3  phil  bll  Ksll  etal  S14 

+  b21  P  Kll  phil  bl2  Ks22  eta2  S24  +  b21  P  Kll  phi2  bl2  Ks21  etal  S14 

+  b21  P  K13  phil  bl2  Ks22  eta2  S24  +  b21  P  K13  phi2  bl2  Ks21  etal  S14 

+  b21  P  Ksl  1  etal  Sll  bl2  Ks22  eta2  S24  +  b21  P  Ksl  1  etal  SI 3  bl2  Ks22  eta2  S24 

+  b22  P  K23  phil  bll  Ksl2  eta2  S24  +  b22  P  K23  phi2  bll  Ksll  etal  S14 
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+  b22  Ks21  eta]  SJl  bll  Ksl2  eta2  S24  +  b22  x?  Ks2I  etal  SB  bll  Ksl2  eta2  S24 
+  b22  Ks22  etal  S21  bll  Ksll  etal  S14  +  b22  Ks22  eta2  S23  bll  Ksll  etal  SmJ) 

>  C2:=  simplify (coeff {collect (expand(zss) ,Md)  ,Md,  1)  )  ; 

C2  :=  {ksII  etal  phi2  bll  Sll  +  Ksll  etal  phi2  bll  SB  +  etal  phi2  i?  SB  bl2  Ks21 

+  etal  phi2  u  SI  1  bl2  Ks21  +  phi2  bll  phi!  u  Kll  +  phi2  bll  phil  u  KB  +  phi2  phil  aB 

+  phi2  phil  bl2  K23  +  phi2  phil  all  +  phi2  phil  bl2  K21 

+  bll  phil  Ksl2  eta2  S21  +  bll  phil  u  Ksl2  eta2  S23  +  phil  bl2  u  Ks22  eta2  S23 

+  phil  bl2  Ks22  eta2  S2l)l{i-bll  Ksll  etal  Sll  b22  Ks22  eta2  S24 

-bll  Kll  phil  b22  Ks22  eta2  S24  -  bll  Ksl2  eta2  S21  b22  Ks21  etal  S14 

-bll  Ksl2  eta2  S23  b22  Ks21  etal  S14  -  all  phil  b21  Ksl2  eta2  S24 

-bll  Ksll  etal  SB  b22  Ks22  eta2  S24  -bll  Kll  phi2  b22  Ks21  etal  S14 

-bll  KB  phil  b22  Ks22  eta2  S24  -  bll  KB  phi2  b22  Ks21  etal  S14 

-  aB  phi2  b21  Ksll  etal  S14  -  all  phi2  b21  Ksll  etal  S14 

-all  phil  b22  Ks22  eta2  S24  -allu^  phi2  b22  Ks21  etal  S14 

-aB  phil  b21  Ksl2  eta2  S24  -  bl2  K23  phi2  b21  Ksll  etal  S14 
-aB  phil  b22  Ks22  eta2  S24  -  aB  phi2  b22  Ks21  etal  S14 

-bl2u^  K21  phil  b21  Ksl2  eta2  S24  -  bl2  K21  phi2  b21  Ksll  etal  S14 

-bl2u^  K23  phil  b21  Ksl2  eta2  S24  +  a21  phi2  bl2  Ks21  etal  S14 

+  a23  phil  bll  Ksl2  eta2  S24  +  a23  phi2  bll  Ksll  etal  S14 

-bl2u^  Ks21  etal  Sll  b21  Ksl2  eta2  S24  -  bl2  Ks21  etal  SB  b21  Ksl2  eta2  S24 

-bl2  Ks22  eta2  S21  b21  Ksll  etal  S14  -  bl2  Ks22  eta2  S23  b21  Ksll  etal  S14 

+  a21  phil  bll  Ksl2  eta2  S24  +  a21  phi2  bll  Ksll  etal  S14 

+  a21  phil  bl2  Ks22  eta2  S24  +  b21  Ksl2  eta2  S21  bl2  Ks21  etal  S14 

+  b21  Ksl2  eta2  S23  bl2  Ks21  etal  S14  +  b22  K21  phil  bll  Ksl2  eta2  S24 

+  b22  K21  phi2  bll  Ksll  etal  SI 4  +  a23  phil  bl2  Ks22  eta2  S24 

+  a23  phi2  bl2  Ks21  etal  S14  +  b21  Kll  phil  bl2  Ks22  eta2  S24 
+  b21  Kll  phi2  bl2  Ks21  etal  S]4  +  b21  KB  phil  bl2  Ks22  eta2  S24 

+  b21  KB  phi2  bl2  Ks21  etal  S14  +  b21  Ksll  etal  Sll  bl2  Ks22  eta2  S24 

+  b21  Ksl  1  etal  SB  bl2  Ks22  eta2  S24  +  b22  i?  K23 phil  bll  Ksl2  eta2  S24 

+  b22  K23  phi2  bl  1  Ksll  etal  S14  +  b22  Ks21  etal  Sll  bll  Ksl2  eta2  S24 

+  b22  Ks21  etal  SB  bll  Ksl2  eta2  S24  +  b22  Ks22  etal  S21  bll  Ksll  etal  S14 
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